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Introduction 


Statement  of  the  problem 

£  <Tl /tJtot  's 

p  Our  motivation  is  to  find  asymptotically  more  accurate  confi¬ 
dence  intervals  for  the  steady  state  mean  of  a  simulated  process. 
By  this  ie  mean^that  the  coverage  probability  error  for  the  con¬ 
fidence  intervals  we  derive  should  be  of  lower  order  than  that  of 
standard  confidence  intervals.  There  are  several  standard  methods 

of  setting  confidence  intervals  in  simulations,  including  the  regener- 

Ht 

ative  method,  batch  means,  and  time  series  methods.  Wb  will  focus^J' 
on  improved  confidence  intervals  for  the  mean  of  an  autoregressive 
process,  and  as  such  our  results  are  useful  outside  of  a  simulation 
setting. 

Improved  methodology  for  setting  simulation  confidence  inter¬ 
vals  is  an  area  of  active  research.  A  recent  survey  article,  Law  AND 
KELTON  (1984),  states 

One  of  the  mosi  important  but  difficult  problems  encoun- 


/  -  T~/-~u .  ■ 
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tered  in  a  real-world  simulation  study  is  that  of  construct¬ 
ing  a  confidence  interval  (c.i.)  for  the  steady  state  mean 
H  of  a  stochastic  proces.  The  information  contained  to 
such  a  c.i.  provides  the  decision  maker  with  a  measure  of 
how  precisely  n  is  known.  Constructing  the  c.i.,  however, 
is  difficult  because  the  output  data  from  a  simulation  are 
in  general  non-stationary  and  autocorrelated,  so  that  di¬ 
rect  application  of  the  techinques  of  classical  statistics  is 
precluded. 

Almost  none  of  the  work  on  confidence  intervals  aims  to  improve 
the  asymptotic  order  of  confidence  interval  accuracy,  however.  If  we 
use  standard  methods  to  set  a  nominal  90  percent  confidence  interval 
for  the  mean  of  the  process,  the  true  proportion  of  the  time  that 
the  actual  mean  falls  within  our  interval  will  be  0.9  plus  an  error 
term  which  is  typically  0(n“l).  The  error  in  half-coverage,  that 
is  the  difference  between  0.45  and  the  true  proportion  of  the  time 
that  the  actual  mean  falls  in  one  half  of  the  nominal  interval,  is  of 
order  0(n-1/a).  Most  work  attempts  to  improve  the  constant  factor 
implicit  in  this  0(n~ l/3).  We  will  derive  a  first  order  correction 
which,  under  reasonable  assumptions,  reduces  the  one  sided  errors 
to  o(n-1/2)  and  a  second  order  correction  which  reduces  these  errors 
to  o(n-1). 

Importance  of  the  problem 

In  the  usual  statistical  context  one  feels  that  a  sample  size  of 
n  =  30  is  reasonably  large — certainly  the  t-distribution  with  29 
degrees  of  freedom  closely  approximates  the  normal  distribution, 
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and  one  regards  normal  theory  as  adequate  for  forming  confidence 
intervals  for  the  mean.  In  forming  a  confidence  interval  for  the 
mean  of  an  asymptotically  stationary  process,  however,  one  may 
need  n  =  10,000  or  greater.  To  illustrate  this,  a  Monte  Carlo  study 
(discussed  in  Chapter  5)  on  the  expected  stationary  waiting  time 
in  an  M/M/1  queue  in  light  traffic  (traffic  intensity  of  0.5)  showed 
among  other  things  that  the  true  probability  that  the  sample  mean 
was  greater  than  the  upper  limit  of  a  90  percent  confidence  interval 
was  not  0.05  but  about  0.14.  Various  criteria  have  been  suggested 
for  evaluating  confidence  intervals,  but  reasonable  accuracy  is  the 
sine  qua  non  of  confidence  intervals. 

For  this  reason,  it  is  desirable  to  find  more  accurate  methods  to 
form  confidence  intervals  for  small  samples.  As  mentioned  above, 
more  accurate  half  coverage  is  desirable.  For  example,  we  might 
wish  to  perform  a  one-sided  test  on  the  estimated  mean. 

Basic  tools  and  assumptions 

We  assume  we  are  given  a  process  x  =  {xf-  :  t  >  1}  which 
satisfies  a  stable  autoregressive  difference  equation, 

(1)  <*o (x<  -M)  +  ••  +  ofc(r,_fc  -  n)  =  e, 

for  i  =  k  +  1, . . . ,  n.  The  e’s  of  the  sequence  {e,-  :  i  >  k  +  1}  are 
zero-mean,  independent  and  identically  distributed  random  vari¬ 
ables  with  moments  of  all  orders  and  a  density  (with  respect  to 
Lebesgue  measure)  which  is  positive  on  an  interval.  We  take  the 
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point  of  view,  then,  that  the  simulator  or  statistician  has  deter¬ 
mined  that  such  a  model  is  appropriate.  If  the  prediction  errors 
of  the  model  (epsilons)  may  only  be  assumed  uncorrelated  and  not 
independent,  the  first  order  correction  still  applies  with  minor  mod¬ 
ifications.  If  the  model  of  actual  interest  is  in  continuous  time,  this 
of  course  necessitates  some  method  like  sampling  at  discrete  inter¬ 
vals.  In  the  regenerative  method  of  forming  confidence  intervals 
for  the  steady  state  value  of  Ef  (x(t)) ,  where  x(<)  is  a  regenerative 
process  in  continuous  time  with  regeneration  times  {U  :  i  >  0},  one 
would  proceed  as  follows.  Let  r,-  be  defined  as  the  regenerative  cycle 
length, 

Ti  =  ti  —  f,-_  i , 

and  let  y,-  be  defined  by 

Sfi=  f  /(*«))<«• 

The  mean  we  want  to  estimate  is  Eyi/Efi,  and  we  may  use  the  Delta 
method  to  derive  the  variance  of  the  natural  estimate  yn/fn »  where 
yn  and  fn  are  the  respective  sample  means.  (See  BILLINGSLEY 
(1979)  for  the  Delta  method).  If  we  are  not  using  the  regenerative 
method,  we  could  sample  the  process  at  fixed  intervals,  which  is 
to  say  that  r,  is  constant,  but  it  may  be  more  convenient  to  sam¬ 
ple  the  process  at  random  intervals.  For  appropriate  choices  of  the 
random  times  {*,•  :  »  >  0}  the  joint  sequence  {(y,-,Tj)  :  t  >  0}  as 
defined  above  will  still  be  asymptotically  stationary  and  will  still 
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obey  the  bivariate  Central  Limit  Theorem,  even  though  it  is  not 
an  independent  sequence.  In  this  case,  Eyt/Efi  is  still  the  natural 
point  estimate,  and  the  Delta  method  may  be  used  exactly  as  be¬ 
fore  to  derive  the  same  variance  constant  as  one  would  derive  using 
the  regenerative  method.  The  only  difference  is  that  the  covariance 
matrix  of  the  process  {(y,-,r,)  :  i  >  0}  is  not  as  easy  to  estimate. 
One  may,  however,  estimate  this  covariance  matrix  using  a  bivariate 
autoregressive  model,  as  in  JOW  (1983),  or  some  other  method.  If  a 
bivariate  autoregressive  method  is  used,  then  the  natural  extensions 
of  the  methods  of  this  thesis  may  be  used  to  obtain  more  accurate 
confidence  intervals,  though  this  extension  to  the  multivariate  au¬ 
toregressive  method  is  not  elaborated  here.  See  FOX  AND  GLYNN 
(1983)  for  further  discussion. 

The  constant  /i  of  (1)  is  the  asymptotic  mean  of  the  x  process, 
and  a0  is  assumed  equal  to  one.  Note  that  the  sequence  of  interest, 
x  =  {ii  :  i  >  1},  need  not  be  stationary,  but  it  will  be  asymptotically 
stationary.  The  requirement  that  the  difference  equation  be  strictly 
stable  is  equivalent  to  the  technical  condition  that  all  roots  of  the 
characteristic  polynomial 


a0zk  +  ••  •  +  akZo 

lie  strictly  inside  the  unit  circle  (see  PRIESTLEY  (1981)). 

If  the  usual  pivot  or  test  statistic  based  on  n  sample  points 
is  tn,  then  the  first  order  pivot  we  will  derive  will  have  the  form 
T\  =  tn  4-  9nn~ l!7  +  pnt\n~l!‘i)  and  the  second  order  pivot  will  be 
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T3  =  Ti  4-  t/ntnn-1  +  unt„n~l.  The  pivots  T\  and  T?  depend  on  n, 
but  we  suppress  this  dependence  in  the  notation.  Ti  differs  from  a 
standard  normal  random  variable  by  an  error  term  which  is  “little 
oh”  of  n-1/3  in  probability,  written  op(n-1/3),  and  7a  by  a  term 
which  is  op(n-1)-  Recall  that  y„  =  op(n~*)  means  that  n* yn  goes 
to  zero  in  probability.  These  test  statistics  will  be  the  basis  of  con¬ 
fidence  intervals  with  coverage  error  probability  of  asymptotically 
lower  order. 

The  basic  tools  behind  our  derivations  are  the  Edgeworth  ex¬ 
pansion,  and  the  Cornish-Fisher  expansion.  KENDALL  AND  STU¬ 
ART  (1977)  contains  a  good  introduction  to  these  concepts.  The 
Edgeworth  expansion  is  an  asymptotic  expansion  for  the  central 
limit  theorem.  If  <J>(z)  denotes  the  standard  normal  distribution 
function  and  <f>(z)  the  standard  normal  density,  the  most  basic  form 
of  Edgeworth  expansion  is 

f’O'-V’EJIw  -Ey)/c  <*}  = 

[-*■  -  *3(*J  -  i)/6-i — ] , 

where  {y,  :  t  >  1}  is  an  independent  and  identically  distributed 
sequence  satisfying  certain  moment  and  smoothness  assumptions. 
The  two  k’s  shown  above  are  0(n-1/3).  Less  basic  forms  of  the 
expansion  will  allow  the  y’s  to  be  non-independent  and  allow  the 
probability  that  the  normalized  sample  mean  lies  in  more  arbitrary 
Borel  sets  to  be  estimated.  In  addition,  we  usually  are  concerned 
not  with  the  normalized  sample  mean  but  with  a  function  of  several 
normalized  sample  means. 
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The  Cornish- Fisher  expansion  amounts  to  a  polynomial  trans¬ 
formation,  which  transforms  the  asymptotically  normal  statistic  tn 
to  a  statistic  with  a  distribution  closer  to  that  of  the  standard  nor¬ 
mal  distribution.  In  fact,  the  first  order  corrected  statistic  T\  above 
is  a  Cornish-Fisher  expansion  of  <n,  and  Tj  is  a  Cornish-Fisher  ex¬ 
pansion  of  Ti . 

Previous  work 

There  are  several  large  areas  of  research  which  provide  a  basis 
for  this  thesis  or  which  are  related  to  it. 

The  first  such  area  is  that  of  Edgeworth  expansions.  Edge¬ 
worth’s  paper  appeared  in  1904.  The  book  BHATTACHARYA  AND 
Rao  (1976)  provides  very  good  background  for  the  independent  and 
identically  distributed  case,  as  does  the  paper  BHATTACHARYA 
and  Ghosh  (1978).  Taniguchi  (1984)  deals  with  the  time  series 
context,  and  GOTZE  AND  HlPP  (1983)  prove  the  validity  of  Edge- 
worth  expansions  for  quite  general  functions  of  weakly  dependent 
random  variables. 

The  original  paper  of  Cornish  and  Fisher  appeared  in  1927.  See 
HlLL  and  Davis  (1968)  for  a  more  recent  study. 

Another  major  area  of  research  is  confidence  interval  methodol¬ 
ogy  in  general.  FISHMAN  (1978)  provides  background  and  presents 
the  autorregressive  method  of  simulation  output  analysis.  Jow 
(1982)  details  the  autoregressive  method  for  vector  processes,  and 
the  articles  Law  AND  KELTON  (1982,  1984)  survey  simulation  ori- 
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ented  work,  and  contain  further  references.  There  is  also  more  sta¬ 
tistically  oriented  research,  such  as  that  of  Efron  cited  below.  The 
electrical  engineering  literature  contains  work  relevant  to  confidence 
intervals,  too.  For  example,  Thomson  (1982)  discusses  sophisti¬ 
cated  methods  for  spectrum  estimation,  and  though  confidence  in¬ 
tervals  are  not  specifically  mentioned,  the  problem  of  estimating  the 
variance  of  the  sample  mean  of  an  asymptotically  stationary  process 
may  be  viewed  as  one  of  spectrum  estimation-see  Jow  (1983)  and 
Priestley  (1981). 

Finally,  there  is  other  work  which  deals  with  improving  the  or¬ 
der  of  confidence  interval  accuracy.  JOHNSON  (1978)  is  perhaps 
the  first  such  article.  He  derives  a  first  order  correction  for  the 
usual  f-statistic  in  an  independent  and  identically  distributed  con¬ 
text.  Glynn  (1982a)  extends  this  to  a  second  order  correction  for 
ratio  estimation,  and  the  methods  of  EFRON  (1984b)  center  around 
the  parametic  bootstrap.  See  also  ABRAMOVITCH  AND  SlNGH 
(1985)  and  the  articles  cited  there. 

Overview 


In  Chapter  2  we  will  develop  machinery  to  help  us  derive  the 
actual  corrections,  though  the  results  of  this  chapter  apply  quite  gen¬ 
erally  to  asymptotically  stationary  sequences.  To  use  the  Cornish- 
Fisher  expansions,  we  need  to  estimate  moments  and  cumulants 
of  some  function  /( i),  where  i  is  a  vector  of  sample  means  of  an 
asymptotically  stationary  sequence.  After  expanding  /  in  a  Taylor 
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series,  we  thus  need  to  estimate  various  joint  moments  and  cumu- 
lants  of  sample  means.  One  aspect  of  this  is  the  need  to  sort  out 
the  contributions  to  a  given  cumulant  in  terms  of  (negative)  powers 
of  the  sample  size,  n.  It  is  primarily  the  second  section  of  Chapter 
2,  “Moment  identities,”  which  details  this  machinery.  The  first  sec¬ 
tion,  “Cumulant  and  moment  bounds,”  facilitates  the  derivations  of 
the  second  section. 

In  this  first  section,  a  result  of  JAMES  (1955,  1958)  and  JAMES 
AND  Mayne  (1962)  is  simplified  and  extended  to  weakly  dependent 
random  variables.  The  result  of  James  asserts  that  the  ;th  cumu¬ 
lant  of  a  polynomial  in  the  sample  average  a0  +  ai(n-1  52?  z*)1  + 

- hafc(n-1  52?  x,)fc  for  independent  and  identically  distributed  x,- 

is  0(nl~3)  as  n  becomes  infinite,  and  this  result  is  referred  to  by 
BHATTACHARYA  AND  GHOSH  (1978)  as  “an  important  combina¬ 
torial  result.” 

Chapter  3  shows  how  to  compute  certain  infinite  sums  associ¬ 
ated  with  sequences  (on  Z+,Z2,  etc.,  where  Z  is  the  set  of  integers 
and  Z+  is  the  set  of  positive  integers)  which  obey  homogeneous  or 
non-homogeneous  difference  equations.  Many  expressions  for  cumu- 
lants  of  functions  of  the  mean  of  an  autoregressive  process  are  of  this 
type. 

Chapter  4  gives  the  actual  derivation  of  the  corrections,  along 
with  an  algorithmic  summary.  Given  the  methods  of  Chapters  2 
and  3  the  derivation  is  relatively  mechanical,  though  by  no  means 
completely  so.  Making  the  Cornish-Fisher  correction  entails  first 
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computing  the  usual  test  statistic,  tn.  One  must  also  estimate  the 
third  and  fourth  moments  of  the  residuals  c,  in  the  autoregressive 
model.  If  n  is  the  sample  size,  estimation  of  covariances  of  the 
underlying  sequence  and  of  the  moments  of  e,-  requires  work  propor¬ 
tional  to  0{n).  The  rest  of  the  computation  requires  work  of  order 
0(k 3),  where  k  is  the  order  of  the  autoregressive  model.  The  first 
order  correction  is  fairly  simple,  and  despite  the  number  of  steps  the 
second  order  correction  is  not  too  computer-intensive. 

Finally,  we  present  numerical  results  in  Chapter  5.  We  consider 
two  genuine  autoregressive  processes  with  independent  and  identi¬ 
cally  distributed  residuals.  In  one  case,  the  residuals  have  a  smooth 
density.  The  second  example  is  the  same,  except  that  the  residuals 
have  a  lattice  distribution.  The  third  model  is  the  M/M/1  queue, 
with  is  not  a  finite  order  autoregressive  process.  The  fourth  model  is 
the  discrete  analog  of  the  M/M/ 1  queue,  that  is,  a  random  walk  on 
the  nonnegative  integers.  This  variety  of  models  is  included  to  test, 
at  least  in  a  few  instances,  the  sufficiency  of  our  sufficient  conditions 
and  the  necessity  of  our  necessary  conditions. 


2 

The  algebra  of 
moments 


In  this  chapter  we  will  develop  tools  which  are  very  useful  for 
simplifying  moment  expressions.  For  example,  to  calculate  the  usual 
test  statistic  (or  “pivot”)  we  need  to  estimate  both  the  mean  of  the 
process  of  interest  and  covariances.  If  the  estimate  of  the  mean  is 
in  and  that  of  the  covariance  at  lag  0  is  /2(0),  in  the  first  order  cor¬ 
rection  we  will  need  to  estimate  the  covariance  of  tn  I?(0).  This 
sort  of  covariance  or  mixed  moment  is  what  we  mean  by  “moment 
expression.” 

First  we  will  obtain  a  result  which  shows,  under  appropriate 
conditions,  that  the  Taylor  expansion  of  a  function  of  a  sample  mean 
has  cumulants  of  the  same  order  as  the  corresponding  cumulants  of 
a  sample  mean.  This  result  and  its  generalizations  will  help  us  to 
derive  the  higher  order  moment  identities  in  the  second  part  of  this 
chapter. 
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Cumulant  and  moment  bounds 

2.1  Introduction.  Let  Sn  denote  the  sum  52*  x»,  where  the  z, 
are  mixing  and  asymptotically  stationary  (in  a  sense  to  be  detailed 
later)  with  asymptotic  mean  zero.  If  we  wish  to  approximate  the 
cumulants  of  v/n(/(5n/n)  —  /( 0)),  we  are  lead  to  consider  the  cumu- 

lantsof  a  Taylor  expansion:  oq  Snn~l/7-\ - hatfcS*nl/2-Ae  =  p(5»). 

In  case  the  sequence  {z,  :  *  >  1}  is  independent  and  identically 
distributed  and  p(5n)  =  Snn~lf 2,  then  the  cumulant  generating 
functions  of  the  polynomial  p  and  of  z  are  related  by  the  formula 

=  ntf)x($n~1/7),  which  shows  that  the  jth  cumulant  for  the 
normalized  polynomial  p,  KJ<n,  is  0(nl~J/7).  JAMES  (1955,  1958) 
and  James  AND  MaYNE  (1962)  have  shown  using  the  machinery 
of  Fisher’s  k-statistics  that  the  same  result  holds  for  general  p(5n) 
when  the  the  sequence  {z,-  :  i  >  1}  is  independent  and  identically 
distributed. 

Our  purpose  in  this  section  is  to  give  a  simple  proof  of  this  fact 
and  some  of  its  extensions. 

Instead  of  considering  the  polynomial  p,  we  take  y/np,  that  is 

<*iSn  H - +  afc5*n1-*,  and  we  will  show  that  the  cumulants  are 

O(n).  We  will  begin  with  the  case  in  which  the  z’s  are  independent, 
but  we  do  not  require  zero  means  or  even  identical  distributions. 

2.2  Review  of  cumulants.  We  record  here  several  facts  about 
cumulants. 
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(1)  The  joint  cumulant  cum(yi, . . . ,  y&)  is  given  by 

U 

where  v  denotes  a  partition  of  {yi,...,yfc}  into  p  sets,  where 
Hui  is  the  expected  product  of  those  y’s  in  the  tth  partition, 
and  where  the  sum  is  taken  over  all  partitions  (p  =  1,...  ,k). 
For  example,  cum(yi,y3)  =  -fityiya  -  EyiEy 3. 

(2)  By  definition,  cum(yi, . . . , y*)  is  the  coefficient  of  •  •  •  a*  in 

the  Taylor  expansion  of  logi?exp[t(3iyi  H - a^y*))  =  Vv(8) 

The  function  0  is  the  cumulant  generating  function.  The  ;th 
cumulant  of  the  univariate  distribution  of  y  is  the  jth  derivative 
the  cumulant  generating  function  at  zero,  D3if/y( 0). 

(3)  The  jth  cumulant  of  y,  which  is  denoted  by  or  */(y),  is  equal 

to  cum(y, . . .  ,y)  (j  times).  This  is  evident  from  the  fact  that 
Kj  =  D3'i>  w  (0)  =  Dl . 3j>r(0)  =  cum(y, . . .  ,y). 

(4)  The  cumulant  cum(xyj,ya, . . . ,  yn)  =  0  if  x  is  independent  of 
the  y’s  and  Ex  =  0,  because  Ex  then  factors  out  of  all  the 
moments  in  the  expansion  of  (1)  above. 

(5)  cum (yi,...,yy)  =  0  whenever  the  y’s  may  be  partitioned  into 
two  groups  independent  of  each  other.  In  this  case,  supposing 
(yi,...,yj)  is  independent  of  (yi+i, . . . ,yy), 

log  2?exp(tiry)  =  log£exp(3iy!  +--+3iyj) 

+  log£exp(j,+  1y,+  1  +  ••  •  +  3yyy), 

whose  Taylor  expansion  has  no  «!  •  •  •  ay  term. 
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(6)  The  ;th  cumulant  of  anj/i  -I - +  aiyi  =  aTy  is  a  linear  com¬ 

bination  of  jth  order  cumulants  which  we  could  write  as 

53  /?i/Cum(y„), 

\v\-i 

where  v  is  a  multiindex,  v  —  |t/|  =  j  means  that 

v  has  j  ordered  elements,  and  cum(y„)  =  c\im(yVl, . . . ,  yVj).  If 
tft i  is  the  cumulac*  generating  function  of  aTy  and  that  of 
y,  then  ^i(a)  =  ^3(are).  Taking  yth  derivatives  on  both  sides 
with  respect  to  a  now  proves  the  assertion. 

(7)  The  cumulant  is  symmetric  and  multilinear.  For  example,  if 
S  =  Sx  Xii 

cum(5,52/n,53/n3)  =  n“3cum(5,52,53) 

=  n-3  53  cum(x,-, ,  x,-,  x,-t ,  x,4  xu  x,-, ) . 

We  are  now  ready  to  prove  the  theorem  for  the  independent 

case. 

2.3  Theorem.  Let  {x,- : »  >  1}  be  an  independent  sequence  and 
let  j  and  k  be  positive  integers.  Suppose  that  the  set  of  moments 
Exi,  •  •  •  x„  is  bounded  for  l  <  jk.  Then  the  j th  cumulant  of  QiS  + 

- 1-  akSnl~k,  where  S  =  x,-,  is  O(n)  as  n  — »  oo. 

Proof.  From  fact  (6)  above,  the  ;th  cumulant  of  ati5  +  •  ••  + 
arj feSnl_fc  is  a  linear  combination  H|i/|=y &/Cum(y„),  where  y  = 
(5, . . . ,  Sknl~k)  and  where  the  coefficients/?,,  do  not  depend  on  n.  It 
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therefore  suffices  to  show  that  each  cum(y„)  is  O(n).  The  expression 
cum(y*,)  means  cum(S‘'1n1-‘'*, . . . , S‘'>nl~v>)  if  u  =  (1/1,.. . 
Notation  may  obscure  the  simple  idea  here,  so  let  us  suppose  we  are 
dealing  with  cum (S,53/n,S3/n,S3/na)  .  That  the  same  argument 
applies  in  general  will  be  obvious. 

Fact  (7)  of  the  previous  subsection  shows  that  this  cumulant  is 

n”4  cum(  xix  *,-4  **, ,  x,-,  xir  ) . 

The  indexes  (*i , . . . ,  »a)  come  from  the  set  of  eight-tuples  of  positive 
integers  less  than  or  equal  to  n,  but  fact  (5)  shows  that  this  cumulant 
is  often  zero.  In  fact,  if  the  four  arguments  split  into  two  indepen¬ 
dent  groups,  the  cumulant  is  zero.  Therefore  we  must  specify  three 
constraints  (the  number  of  arguments  minus  one)  to  make  it  non¬ 
zero.  That  is,  if  each  argument  is  a  node  and  if  we  regard  two  nodes 
as  connected  whenever  the  arguments  are  dependent,  we  must  have  a 
connected  graph.  Two  arguments  in  our  context  are  dependent  only 
when  they  share  an  index — for  example  zn  and  z„Xj,  are  dependent 
only  if  =  »a  or  %i  =  t3.  Note  that  the  possible  number  of  such  con¬ 
straints  does  not  depend  on  n.  Thus  the  cumulant  is  nonzero  only 
if  the  indexes  satisfy  at  least  one  set  of  three  equality  constraints 
from  a  fixed  number  of  j  possible  sets.  It  follows  that  the  number 
of  nonzero  terms  in  the  sum  J^cu m(z,-,,ZiJz,t,z,-4zlt,Ziiz,Tz,i)  is 
0(n8-3)  =  0(n5).  Because  of  our  boundedness  assumption  on  mo¬ 
ments,  and  taking  into  account  the  denominators  in  the  cumulants’ 
arguments,  the  total  sum  is  0(n5/n4)  =  O(n)  as  required. 
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In  general,  each  new  block  x3l  ■  •  •  x3 ,n1-<  which  is  an  argument 
of  the  cumulant  increases  the  size  of  the  space  of  possible  indexes 
by  a  factor  of  n! ,  but  by  introducing  one  additional  necessary  con¬ 
straint  and  dividing  by  n*-1,  the  growth  of  the  sum  remains  at 
O(n).  Furthermore,  if  there  is  only  one  block,  no  constraints  are 
needed.  □ 

Before  proving  the  theorem  for  non-independent  sequences,  we 
will  review  the  concept  of  mixing. 

2.4  Review  of  mixing.  The  useful  hypothesis  of  mixing  gen¬ 
eralizes  the  idea  of  independence.  Let  {xt-  :  t  >  1}  be  a  sequence 
of  random  variables  and  suppose  the  event  A  is  in  the  sigma  field 
<t(x  i,...,ifc),  while  B  €  <r(xfc+», . . . ).  Then  if  the  absolute  differ¬ 
ence  of  probabilities|Pr(>4B)  -Pr(i4)Pr(JB)|  is  less  than  or  equal  to 
an  uniformly  for  all  such  A  and  B,  we  say  that  the  x’s  are  mixing, 
provided  an  — »  0.  The  an  are  referred  to  as  mixing  constants. 

Even  to  prove  a  central  limit  theorem  for  the  x’s,  we  would 
require  mixing  or  something  which  implies  it.  See  for  example 
Billingsley  (1968),  pages  166  and  174.  One  may  also  show 
(BILLINGSLEY  (1979),  p.  317)  that  if  the  random  variable  y  is 
measurable  <r(ii , . . . ,  x*)  while  z  is  measurable  a(ifc+n, . . . ),  and 
if  y  and  z  have  bounded  fourth  moments,  then  \Eyz  —  EyEz\  < 
8(1  -t-  Ey*  +  Ez*)aJ2 . 

2.5  Theorem.  Let  {x,-  :  »'  >  1}  be  a  sequence  such  that  (1) 
for  some  positive  integers  j  and  k,  the  set  of  moments  Eiil  ■  ■  x,-( 
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is  bounded  for  l  <  4(jk  —  1)  and  (2)  the  sequence  is  mixing  with 
mixing  constants  an  =  0(n~2^3~l+e))  for  some  e  >  0.  Tien  the 

jtb  cumulant  of  a-i 5  + - h  afcSn1-fe,  where  S  =  2»>  JS  0(») 

as  n  — »  oo. 

Remark.  In  applying  the  lemma  of  BILLINGSLEY  cited  in  the  last 
subsection,  the  random  variables  y  and  z  will  be  products  of  at 
most  kj  -  1  z’s,  say  y  =  x„  —  xUl  •  ••  x„m ,  and  similarly  z  =  xp. 
Then  the  lemma  together  with  the  hypotheses  of  our  theorem  guar¬ 
antee  \Exvxp  —  ExvExp\  is  uniformly  0(g~^3~l+(^)  whenever  the 
minimum  p,-  exceeds  the  maximum  i/f-  by  at  least  g  . 

Proof.  As  in  Theorem  2.3  we  only  need  to  show  that  each  cum(y„) 
is  O(n),  and  this  reduces  to  considering  the  same  sort  of  sum  of 
cumulants.  The  only  difference  is  that  terms  like 

y]  cum(xtlJ  i| ,  *,•, ,  z,4  z,# ,  x,-,  x,T  x,-g ) 

may  never  be  zero,  but  with  a  strong  enough  mixing  condition  the 
sum  of  such  terms  will  still  be  0(n),  as  we  shall  see. 

To  be  specific,  suppose  again  that  we  are  dealing  with 

n  y  1  cum(xi4 ,  X|S  x,j,  zwx,-5 ,  x,-,x,Tx,-t ) 

which  we  want  to  be  0(n).  We  have  t  =  4  arguments  of  the  cu¬ 
mulant,  d  =  8  dimensions  to  the  space  of  indexes,  c  =  £  —  1  =  3 
necessary  constraints  for  a  nonzero  cumulant  in  the  independent 
case,  and  p  =  4  total  powers  of  n  in  the  denominator,  which  have 
been  factored  out  of  the  sum.  The  relation  d  —  c  -  p  =  1  always 
holds,  and  this  guarantees  that  the  sum  is  0(n)  in  Theorem  2.3. 
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Suppose  that  <■••<*«  =  id-  It  is  enough  to  show  that  the 
sum  over  such  terms  is  O(n),  as  long  as  we  can  do  the  same  for  each 
ordering  of  the  t’s.  We  will  need  to  consider  the  cth  smallest  gap 
%k  -  ik-i  in  the  index  »  =  (*!,...  ,ia),  which  we  call  g.  If  the  gap  g 
is  “large",  then  x’s  whose  indexes  differ  by  g  or  more  will  be  nearly 
independent.  If  such  x’s  actually  were  independent,  and  if  the  cth 
smallest  gap  in  (t'i , ... ,  »<j)  is  g ,  then  the  corresponding  cumulant 
would  be  zero.  In  our  specific  example  (c  =  3),  this  means  that 
only  two  pairs  of  x’s  or  one  triple  could  be  dependent,  which  in  turn 
implies  that  one  of  the  cumulant ’s  arguments  would  be  independent 
of  the  rest,  exactly  as  in  the  proof  of  Theorem  2.3.  With  our  mixing 
condition,  the  cumulant  will  be  nearly  zero.  Therefore,  we  must 
answer  two  questions: 

(1)  If  the  cth  smallest  gap  in  (tj, . . . ,  »<*)  is  g ,  then  how  small  is  the 

corresponding  cumulant,  and 

(2)  How  many  points  (t'i ,...,*<*)  have  cth  smallest  gap  equal  to  g ? 

To  answer  the  first  question,  we  expand  the  cumulant  in  terms 
of  moments,  as  in  (2.2.1).  If  we  factor  each  expected  product  in  the 
moment  expansion  as  though  the  x’s  whose  indexes  differ  by  g  or 
more  were  independent,  we  arrive  at  an  expression  which  is  identi¬ 
cally  zero  (because  g  is  the  cth  smallest  gap).  But  now  each  time  we 
factor  such  an  expected  product  we  incur  an  error  which  is  uniformly 
0(<r<>-l+<>)  — see  the  remark  preceding  the  proof  of  this  theorem. 
Note  that  each  term  in  the  cumulant  expansion  is  the  product  of  at 
most  d  moments,  and  each  moment  is  the  expected  product  of  at 
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most  d  x’s.  If  b  >  I  is  our  uniform  bound  on  x-moments  of  order 
d  or  less,  then  each  term  of  the  cumulant  expansion  is  bounded  by 
bd  in  absolute  value.  There  are  exactly  d  x’s  in  each  term  of  the 
expansion  (each  index  occurs  once),  and  each  time  we  factor  we  in¬ 
duce  an  error  bounded  by  0(g~^3~l+e^).  The  total  absolute  error 
induced  by  factoring  expected  products  is  therefore  0(g~^3~l+t^). 
This  answers  the  first  question:  cumulants  whose  cth  smallest  gap  is 
g  are  uniformly  0(g~^3~1+f^)  in  absolute  value  ( before  multiplying 
by  the  factor  n~p  =  n~*  in  our  example). 

Next  we  must  bound  the  number  of  t<j)  with  cth  smallest 

gap  equal  to  g.  Recall  that  ij  <  •  •  •  <  i8  =  ij,  and  let  gi  = 
*i, . . . , gd  =  id  —  t'd-i-  The  cth  smallest  gap  is  the  cth  smallest  of 
<fa, . . .  ,g<i-  Again,  it  suffices  to  find  a  bound  based  on  the  assumption 
that  <fa, . .  • ,  ge+x  are  less  than  or  equal  to  yc+3,.  •  •  ,9d-  There  are 
g c  ways  of  assigning  values  in  <7}  to  the  c  smallest  g’s  and 

(g  —  l)c  ways  of  assigning  the  values  {1,. . .  ,g  —  1}.  Therefore  there 
are  0(gc~l)  ways  to  assign  g}  such  that  the  cth  smallest 

gi  is  exactly  equal  to  g.  Some  of  these  ways  may  not  satisfy  the 
constraint  £2  9i  —  but  that  does  not  affect  our  upper  bound. 
There  are  fewer  than  nd~c  ways  to  assign  the  remaining  <7,-’s,  hence 
the  number  of  (tj, . . . ,  14)  with  cth  smallest  gap  of  g  is  0(gc~l nd~c). 
This  answers  the  second  question  stated  on  above. 

We  know  already  that  the  number  of  cumulants  in  the  sum  with 
g  =  0  is  exactly  nd~e.  Therefore  the  entire  cumulant  sum,  after  di¬ 
viding  by  the  denominator  np  =  nd~c~l  is  bounded  in  absolute  value 
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by  a  constant  times  n~(d~e~l\nd~e  4-  nd~c  9e~1  /93~1+e)  = 

0(n)  if  j  —  1  +  e  >  c.  For  ;th  order  cumlulants,  the  number  of 
constraints  c  is  j  —  1,  so  the  cumulant  sum  is  0(n).  □ 

2.0  Remark.  Note  in  particular  that  the  ;th  cumulant  of  S, 
cum(5, ...,5)  ,  is  O(n),  where  S  =  x*>  provided  the  mixed  x- 

moments  of  order  4(j  —  1)  are  uniformly  bounded  (we  assume;  >  1). 

We  also  have  cum(5, . . . ,  S,  xn)  =  0(1),  since  the  only  con¬ 
tributing  terms  are  those  for  which  the  indexes  ti,...,ty  are  all 
approximately  equal  to  n.  In  general, 

cum(5, . . . ,  5,  *»-.,)  =  0(1). 

Here  the  it  are  fixed  nonnegative  integers,  not  necessarily  distinct. 
For  this  to  hold  (assuming  5  repeats  j  times)  it  suffices  that  the 
mixed  x-moments  of  order  4(j  +  m  —  1)  be  uniformly  bounded  and 
the  mixing  constants  am  be  of  order  0(m'"2*JI+<l). 

In  this  connection  we  have  the  following  corollary. 

2.7  Corollary.  Let  S  =  x,.  Then  for  positive  integer  k  the 

moment  ESkxn+il  ■  ■  ■  xn+,r  is  of  order  np,  where  p  is  at  most  [k/ 2], 
provided 

(1)  E\xi\  =  0(n"2). 

(2)  The  kth  cumulant  of  S  is  0(n). 

(3)  cum(S,.. .  S,xn_j,  •  •  •  xn-»,)  =  0(1)  (S  repeats  A:  or /ewer  times 
as  an  argument). 

(4)  x  has  uniformly  bounded  mixed  l  +  k-order  moments. 
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Proof.  Use  induction  on  k  +  l.  If  k  +  l  <  1,  the  conclusion  follows 
immediately  from  (1)  and  (4). 

For  Jfc  >  2  and  l  =  0  we  observe  that  the  Ath  cumulant  of  a 
is  0(n),  expand  the  cumulant  in  terms  of  moments,  and  apply  the 
induction  hypothesis.  The  first  term  in  the  cumulant  expansion 
is  ESk,  while  the  other  terms  are  products  \[ESVi  with  = 
A — so  the  induction  hypothesis  proves  this  case.  For  general  k  + 
l  >  2,  consider  cum(5,...,5, xn_,t  ,  which  is  0(1)  by 

hypothesis.  □ 

2.8  Remark.  In  a  sense  these  results  may  not  seem  as  good 
as  one  might  hope,  and  they  can  in  fact  be  strengthened.  Sup¬ 
pose  we  are  interested  in  the  second  cumulant  of  S/n  +  S3/n3. 
Then  we  have  shown  that  cum(S/n,  S/n)  ,  cum(S/n,S3/n3)  .  and 
cum(52/n3,53/n3)  are  all  0(l/n).  This  means  that  the  least  sig¬ 
nificant  term  in  the  Taylor  expansion  has  as  much  impact  as  the 
most  significant  term.  We  can  sharpen  the  result  if  the  sequence 
{x,-  :  i  >  1}  has  asymptotic  mean  zero,  which  will  normally  be  true 
in  the  cases  of  interest  to  us. 

2.0  Corollary.  Let  {x,  :  i  >  1}  be  an  independent  sequence  of 
random  variables  with  zero  mean  and  uniformly  bounded  moments 
of  order  p  =  p,-  (the  p<  are  positive  integers ),  and  let  S  =  x,-. 

Then  cum(5Pl, . . .  ,SP>)  is  0(nlp/3I)  and  also  0(np~J+l).  ftp/ 2] 
means  the  greatest  integer  in  p/2.) 

Proof.  The  0(np~,+l )  bound  follows  from  the  argument  of  Theorem 
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2.3:  there  must  be  j  —  1  constraints  on  the  indexes  of 

cum (llj=  i  -  ) 


for  a  nonzero  cumulant.  It  is  also  true,  however,  that  each  index  »m 
must  be  equal  to  at  least  one  other  index.  If  for  example  the  index 
im  appears  only  once,  then  Ex<m  =  0  factors  out  of  each  term  in 
the  expansion  of  the  cumulant  as  a  sum  of  moments.  Thus  for  a 
nonzero  contribution  we  require  \p/ 2]  constraints,  and  the  sum  is 
therefore  0(nlp/2l).  □ 

Either  of  the  two  bounds  in  the  statement  of  the  corollary  may 
be  sharper.  Our  new  bound  shows,  for  independent  zero  mean  z’s, 
that  cum (S/n,52/n2)  and  cum(52/n2, 52/n2)  are  0(n-2),  which 
is  an  improvement. 

The  simplest  asymptotically  stationary  case  occurs  when  the  z’s 
are  independent  with  zero  mean,  except  that  Exi  ^  0.  In  this  case 
for  a  nonzero  cumulant  each  index  of  cum  (n?=i  *.»»•  ••)  must 
be  equal  to  another  index  or  equal  to  one,  and  this  still  requires 
[p/2l  constraints.  In  the  dependent  asymptotically  stationary  (zero 
mean)  case,  this  means  -!ach  index  must  be  near  another  index  or 
near  zero. 

For  how  many  sets  of  indexes  (*i , . . . ,  *p)  is  each  index  between 
one  and  n  and  also  within  g  units  of  zero  or  of  another  index?  We 
can  bound  this  number  by  cnlp/2l< f  where  c  is  independent  of  n 
and  g.  Therefore  the  number  with  this  minimum  gap  equal  to  g  is 
bounded  by  c,nlp/2l^p_l.  Suppose  2£|xm|  =  0(m~l)  and  the  mixing 
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constants  am  are  0(m  2*).  We  then  find  that  the  cumulant  of 
the  corresponding  selection  of  z's,  namely  cum^ffyii  *s 

0(g~l). 

Combining  these  observations  with  the  argument  of  Theorem 
2.5  gives  the  following. 

2.10  Corollary.  Let  p  =  p,-,  where  each  p,  is  a  positive 

integer,  and  S  =  x,  (p  >  2).  Suppose  that 

(1)  mixed  moments  of  x  *s  of  order  4 (p  —  I)  are  uniformly  bounded. 

(2)  The  x’s  are  mixing  with  mixing  constants 

am  =  0(m-2<p+<>)  (e  >  0) 

(3)  £|xm|  =  0(m-<P+0). 

Then  cum( 5Pl , . . . , )  is  0(n\*l7\).  □ 

Note  that  the  argument  of  Theorem  2.5  implies  that  the  cumu¬ 
lant  above  is  0(np_,+  l)  if  conditions  (2)  and  (3)  above  are  replaced 
by  the  requirement  that  the  mixing  constants  om  =  0(m_2^~1+el). 
We  need  to  make  two  important  remarks. 

2.11  Remark.  Suppose  some  or  all  of  the  arguments  to  the 

cumulant  are  multiplied  by  a  product  of  rn_f’s.  In  other  words, 
each  argument  is  Sp>  Il/s’i  where  we  allow  the  product  to 

be  empty  or  py  =  0  (but  not  both).  If 

(1)  The  sum  £]py  is  denoted  by  p, 

(2)  the  total  number  of  indexes  is  d  =  p  +  J2  mj. 

(3)  the  number  of  arguments  for  which  mj  =  0  is  j, 
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(4)  at  least  one  my  >  0, 

(5)  the  i-sequence  has  bounded  mixed  moments  of  order  4(d  —  1), 
and 

(6)  the  mixing  constants  am  =  0(m-2^+<l)  (e  >  0) 
then  the  cumulant  is  0(np~3). 

This  is  an  extension  of  Theorem  2.5.  The  corollary  above  can 
be  similarly  generalized,  but  we  shall  not  do  so. 

2.12  Remark.  Given  a  finite  number  /  of  asymptotically  sta¬ 
tionary  sequences,  a:*1), . . . ,  x W  ,  we  may  take  =  Yassi  x\3^  and 
obtain  all  the  analogous  results  for  mixed  cumulants  provided  the 
vector  process  {(zjl\ . . . ,  zj^ )  :  *  >  0}  satisfies  the  mixing  and  mo¬ 
ment  conditions. 

We  will  come  across  some  estimates  of  the  form 


Here  we  are  simplifying  slightly,  because  the  sums  5/n  and  the 
factors  z*/n  could  be  from  different  sequences,  and  furthermore  we 
might  actually  have  products  of  factors  of  the  form  x, *-»y/n,  with 
a  different  j  for  each  factor.  A  primary  example  is  the  estimate  of 
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covariance.  We  have 

1  n_1 

Rl  =~  ( xi  ~  %n)(xi+l  —  ®n) 

”  1 

1  n"1 

=  “  2 (*•'  “/*)(*«'+ 1  “/*) 

U)  *,  w  » 

“  -(*»  “/*)(*»+!  -M) 

Tl 

-  (*.  -  (i + i) 

+  m)  [(r, 

There  are  two  sequences  involved  here:  (*, -m)  and  (*,■  -/i)(xt+i  - 
/x) .  Each  term  does  have  the  form  suggested  above,  except  that 
there  may  be  an  extra  factor  of  n_p.  We  want  to  conclude  that  the 
cumulants  of  such  expressions  have  the  same  asymptotic  behavior 
as  do  the  Taylor  expansions  dealt  with  in  Theorem  2.5. 

Theorem  2.5,  after  rescaling,  shows  that  has 

cumulants  «y  of  order  0(n1-J).  This  ;th  cumulant  is  a  linear  com¬ 
bination  of  terms 

cum(5Pl, . . . ,  Sp>). 

Subject  to  appropriat;  mixing  and  moment  conditions,  Remark  2.11 
shows  that  if  an  (xn/n)g  argument  is  adjoined,  then  the  order  of  the 
cumulant  is  actually  reduced.  Alternatively,  if  a  term  Sp>  is  multi¬ 
plied  by  (zb/h)*  for  integer  q  >  0,  then  the  order  of  the  cumulant  is 
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unchanged  if  q  =  1,  and  otherwise  reduced.  Therefore  we  have  the 
following  theorem. 

2.13  Theorem.  Let  /„  =  X^_0(5/n)m(zn/n)«-"crm,  where  the 

qm  are  nonnegative  integers  and  S  =  Yl’i  x*‘  p  =  j  •  [maxm(m  + 

gm)]  >  1.  Assume 

(1)  the  z’s  have  uniformly  bounded  mixed  moments  of  order  4 (p  - 
1),  and 

(2)  the  x’s  are  mixing  with  mixing  constants  am  = 

Then  the  j th  cumulant  of  /„  is  0(nl~J)  as  n  — »  oo.  □ 

2.14  Remark.  Suppose  a  finite  number  of  sequences  x^ , . . . ,  x ^ 
satisfy  the  vector  versions  of  (1)  and  (2)  in  the  last  theorem.  Here 
(1)  means  that  the  moments  of  form  Exj™1^  ■  ■  ■  are  uniformly 
bounded  for  l  <  4(p  —  1),  and  (2)  becomes  the  obvious  extension  of 
the  mixing  concept.  We  may  then  replace  any  factors  (5/n)  by  any 

n-1 1C"  and  any  xn/n  by  zjj^/n  for  fixed  i,  or  by  x\ m^/n  for 
fixed  i.  Also,  dividing  some  terms  by  additional  factors  of  n  certainly 
will  not  increase  the  order. 

Let  us  illustrate  how  the  estimate  Ri  of  equation  (2.12.1)  fits 
this  form.  Let  x ^  be  the  sequence  (z,-ji)  and  z ^  be  the  sequence 
(xi  -  p)(z<+l  -  p).  Then  the  standard  estimate  of  (2.12.1)  is 


SW 

n 


We  can  therefore  conclude  that  the  jth  cumulant  of  Ri  is 
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0(n1_'>).  The  maximum  number  of  x  factors  in  any  term  of  R\ 
is  two,  so  it  will  suffice  to  have  bounded  mixed  x  moments  of  order 
4(2;  -  1)  and  mixing  constants  am  =  0(m-3*J+e*). 


Moment  identities 

2.15  Conditions.  The  identities  in  this  section  require  several 
types  of  conditions  on  the  sequence  {xi  :  i>  1}.  The  parameters  p 
and  l  below  are  specified  in  each  subsection. 

(1)  Mixing.  The  z’s  should  have  mixing  constants  an  =  0(n-3p). 

(2)  Moment  condition.  Mixed  x-moments  E |xtl---xt>|  are  uni¬ 
formly  bounded  for  all  integer  ;  less  than  or  equal  to  the  order 
of  the  moment  being  considered. 

(3)  Asymptotic  stationarity.  If  E  denotes  expectation  and  E,  de¬ 
notes  expectation  with  respect  to  the  limiting  stationary  distri¬ 
bution  (assumed  to  exist),  then  | E,Xix  •  •  •  x,m  —  2?x,-,  •  •  •  x,m|  = 
0(n~p)  whenever  the  smallest  index  ij  is  at  least  n  and  the 
number  of  indexes  (m)  is  no  more  than  l. 

(4)  The  asymptotic  mean  is  zero. 

2.1®  Notation.  We  will  be  expressing  quantities  like  E(%2*  *»)p 
as  polynomials  in  n,  and  we  need  to  establish  a  notation  for  the 
coefficients.  “A”  in  the  notation  below  is  meant  to  indicate  that  the 
asymptotic  means  have  been  subtracted  off,  so  in  effect  the  x’s  are 
assumed  to  have  asymptotic  mean  zero.  If  Sx  denotes  J2*  x«»  we 
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will  write,  for  example, 


=  p0(ASx)nQ  +  O(f(n)) 


=  p  1(A5,a)n1  4-  M°(A5s2)n°  +  0(/(n)) 
=  M1(A5,3A5v)n1  +  M°(AS,2AS,)n° 


+  0(/(n)) 


£(fly  -  fly)4  =  [/i-2(Afly4)n-2  +  M~3(Afly4)n-3 
+  M-4(Afly4)n-4  +  0(/(n))] 


and  so  on.  The  0(/(n))  above  is  a  generic  error  term  which  will  be 
specified  in  what  follows. 

Thus  the  superscript  of  p  corresponds  to  the  power  of  n,  and 
the  argument  of  pp(-)  is  the  quantity  whose  moments  are  being 
considered.  The  notation,  then,  decomposes  the  expectation  into 
its  components. 

We  will  write  p,  to  refer  to  expectations  with  respect  to  the 
limiting  stationary  distribution.  In  addition,  if  fl  is  a  vector  of 
estimates  and  a  is  a  vector,  pp(ARTa)  is  2ZMp(Afl,o»),  and  a 
similar  notation  applies  to  matrix  estimates  M. 

If  the  context  is  clear,  we  may  drop  the  argument:  pp.  The  hat¬ 
ted  notation  pp{  )  refers  to  the  estimate  of  the  un-hatted  quantity. 
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We  may  also  need  to  consider 


E 


and 


E 


Fortunately,  in  such  contexts  the  underlying  sequence  x <  will  be 
evident,  and  we  will  use  the  same  notation  as  before  but  with  argu¬ 
ments  (A5f  AS'I)  and  (ASP Ax*+1),  respectively. 

2.17  First  moments.  Evidently, 


(1) 


»5(AS.)  =0 

oo 

M°(A  S,)=X>x,-, 


with  error  term  of  order  0(n  (p  *1),  provided  p  >  1  in  the  condi¬ 
tions  of  Subsection  2.15. 

2.18  Second  moments.  Given  that  Sn  =  x,-,  and  Ri  is  the 

correlation  of  the  z  sequence  at  lag  »,  in  the  stationary  case  we  get 


-ESi 


) 
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«—  1  f  00  00  \ 

=  12  (Jt°  +  12  ^  ~ 2  12  ^•) 

k=0  '  1  fe-J-1  ' 

(OO  v  00  00 

5o  +  2j^^,j-2^t5,+2  £  (t  -  n)i?,, 

1  '  i  «=»+i 

provided  .S*  is  0(n-p)  for  p  >  2,  in  which  case  the  remainder  term 
(the  last  term  above)  is  0(n-lp-2l). 

The  asymptotically  stationary  case  is  similar,  but  there  is  an 


additional  term 


^  ^  Et i.ifc. 


k=l t=l 


The  subscripted  E,  denotes  expectation  with  respect  to  the  asymp¬ 
totically  stationary  distribution. 

From  this, 


M1(A5„)=Mi(A5n)  =  /Z0  +  2^^ 

1 

OO 

rf(AS„)  =  -2£.R, 


P°(AS„  )  =  /!«( AS»)  +  £  £  **<*»  -  *.*<*»• 


fc=l 1=1 


To  estimate  the  error  term  in  the  above,  we  need  to  bound  the 


tail  of  the  double  sequence,  which  is 


12  12  ExiXi  ~  E»x'xi  ~1212  Exixj  -  E, 
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Applying  these  techniques  leads  to 

M1(A5n3)=  £  EtxoXiXj  +  3/i°(A5„)/iJ  (AS*2). 

t,j=-oo 

The  conditions  are  that  p  >  2,  and  l  >  8  in  the  conditions  of  Sub¬ 
section  2.15.  We  will  not  need  /i°(A5n3). 

2.20  Fourth  moments.  We  now  can  begin  to  use  the  cumulant 
bounds  developed  in  the  first  part  of  this  chapter.  If  xn  satisfies 
the  hypotheses  of  Theorem  2.5  for  fourth  order  cumulants,  we  can 
expand  the  fourth  cumulant  of  5  =  xi  *n  terms  of  moments: 

O(n)  =  cum(5,  S,S,S) 

=  ES 4  -  4ES*ES  -  SETS7  +  liES'&S  -  6E*S. 
Equating  the  0(n3)  terms  of  the  above  shows  that 

In  view  of  the  more  general  results  of  the  first  part  of  this  chapter, 
we  can  certainly  draw  a  similar  conclusion  for  quantities  like  Rj.  We 
will  also  need  to  know  ^l(A$n4).  In  the  context  of  autoregressive 
processes,  we  will  have  an  easy  way  of  estimating  the  coefficient  in 
the  stationary  case,  but  we  will  need  the  correction  for  the  asymp¬ 
totically  stationary  case.  The  method  of  derivation  is  the  same  as 
that  indicated  in  the  last  subsection.  After  some  algebra  we  find 

M‘(ASn4)  =  -3  (n\  (A5n3))2  4-  E,x*  -I-  4^  (AS  Ax3+1) 
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+  6m?(A52Az2+1)  +  4M2(AS3A*i+1) 

+  V(A5n)/*i(AS„3) 

+  6  [m°(A Sn2)  -  m2(A5„2)]  (AS„2). 

The  last  two  lines  above  gives  the  correction  for  the  asymptotically 
stationary  case.  For  the  coefficients  of  this  subsection  to  be  valid,  it 
suffices  that  l  >  12  and  p  >  4  in  Subsection  2.15. 

2.21  Fifth  and  sixth  order  moments.  If  xn  satisfies  the 
hypotheses  of  Theorem  2.5  for  fifth  order  cumulants,  then 

M2(AS„5)  =  10M1(A5n3)/iI(^^n2)  -  15  (nl(ASn))2  n°(ASn) 
M3(A5„fl)  =  15  (/x1  (A5n2))3  . 

It  is  fortunate  that  not  only  are  the  most  significant  terms  of  the 
higher  order  moments  reduced  to  coefficients  of  lower  order  mo¬ 
ments,  but  also  that  the  nonstationary  correction  is  based  only  on 
first  and  second  moment  nonstationary  terms. 

In  the  next  subsection,  we  extend  these  formulas  to  the  case 
of  mixed  moments.  These  results  will  be  very  useful  in  subsequent 
calculations. 

2.22  Mixed  moments.  Given  two  or  more  sequences,  {a,  : 
t  >  1},  {6,-  :  t  >  1},  and  so  on,  we  may  easily  extend  the  moment 
identities  of  the  preceding  paragraphs  by  identifying  the  coefficients 

of 


E{3Y,a<  +  t5I6<)  ’ 
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provided  sa,-  +  tb<  satisfies  the  sufficient  conditions  of  Subsection 
2.15.  This  gives  us  the  formulas  enumerated  below,  where  Sa  = 
53?  a,-,  and  similarly  for  Sb  and  Sc. 

In  many  cases,  a  “/*”  may  be  replaced  by  a  Those  formulas 
derived  using  Theorem  2.5  (which  means  the  fourth  order  and  higher 
ones)  will  apply  to  non-sample  means  like  Rj. 

OO 

(1)  fi'i&Sa&Sb)  =  E$aQbQ  +  ^2(Esaobi  +  £,at60) 

1 

OO 

(2)  n°(i\Sai\Sb)  =  -  ^fc(£,a06fc  +  E.akbo) 


+  ( Eat  bj  —  Efdibj ) 


Hl  (ASa2  ASb)  =  ^2  Eta0aibj 


t,;  =  — oo 


+  2n°(*Sa)nl(ASaASb) 

+  n\  (ASa2)n°(ASb) 

(4)  /x2(A5a2AS6A5c)  =/ii(ASa2)/zJ(A56ASc) 

+  2^,(ASaASb)ti\(ASaASc) 

(5)  M2(A5a4A56)  =  4M1(A5a3y(ASaA56) 

+  6nl(AS<*ASb)nl  (A  Sa2) 

-  3  (a*1  (A503))3  /i°( A5(,) 

-  12/il  (ASa2)/!1  (A50  AS(,)/i0(A5a 


Summing  difference 
recursions 


In  our  moment  calculations  we  will  have  to  compute  sums  like 

OO  00  oo 

X)y*»  or 

0  0  i=0 

where  the  y’s  satisfy  a  difference  equation  and  the  to’s  satisfy  a 
difference  equation  in  each  index. 

Sequences  indexed  by  non- negative  integers 

3.1  One  dimensional  recursions.  We  assume  that  the  sequence 
{y» :  *  >  0}  satisfies  a  strictly  stable  order  k  difference  equation, 

y»  +  °iy»-i  h — Ofcyn-*  =  o 

for  n  >  k,  where  yoi-  ->y*-i  are  known.  That  the  recursion  is 
strictly  stable  means  that  all  the  roots  of  the  characteristic  polyno¬ 
mial  zk  -I-  +  •  •  •  +  OfcS0  lie  in  the  interior  of  the  unit  circle. 

Under  this  assumption,  the  sums  in  question  converge  absolutely. 
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(See  PRIESTLEY  (1981)  for  further  discussion  and  results).  We 
may  calculate  these  sums  directly,  or  by  using  the  generating  func¬ 
tion  <p{z)  =  Vizi-  The°  =  £o°  y«  and  y>'(l)  =  jy y. 

Because  the  y,’s  satisfy  the  homogeneous  difference  equation, 
we  find 

<p{z)(  l  +  •  •  •  +  afczfc)  =  yo  +  «(yi  +  «»iyo)+ 

— p  *fc_1(yfc-i  h - h  flfc-iyo) 

and 

=  VW 

o 

m  -  yot1  -b--  +afc-i)  +  ••  •  +  yjb-i(i) 

(i +  •••  +  «») 

It  is  possible  to  calculate  £  j'yy  by  considering  <p'{  1),  but  practically 
speaking  there  is  an  easier  method.  We  have 

oo  oo  i 

x  y=i  1=1 

00  oo 

«=i  j=« 

Note  that  the  sums  5,-  =  yy  obey  the  same  difference  recursion 
as  the  y,-’s,  but  with  different  initial  values.  Therefore  a  convenient 
method  to  compute  this  sum  is  as  follows. 

3.2  Algorithm.  To  compute  YLT  3  V v 
(1)  Compute  So  =  V*  fr°m  (I)  above. 
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(2)  Compute  Si  =  S0  -  yo, ■  ■  ■  ,Sk  =  Sk-\  -  yk- 

(3)  The  answer  is  then 


Si  (1  +  •  •  •  -f  afc-i)  +  •••-(-  5fc(l) 
(1  -I - +  a*) 


□ 

3.3  Two  dimensional  recursions.  We  will  also  need  to  sum 
two  dimensional  arrays  whose  elements  satisfy  a  difference  equation. 
Specifically,  if  J2,y  is  the  non-stationary  expectation  E(xi— n)(ij  —  fi) 
and  Ri  is  the  stationary  expectation  Et(ii  —  p)(x0  —  /*),  we  will 
need  to  sum  the  array  Rij  —  R\i-j\  in  its  entirety  and  also  along  its 
diagonals. 

In  an  autoregressive  model  of  order  k  in  which  the  sequence 
{x i-n  :  *  >  1}  satisfies  the  usual  difference  equation  and  y,-  = 
we  get 

Rij  =  Eyiyj 

=  -Eyi(aiyj-i  +  ■■■  +  afcyy_fc  +  €y) 

=  —  (aiRij-i  ^ - 1-  akRitj-k  +  Eyitj) , 

assuming  j  >  k .  This  gives  Rij  as  the  solution  of  a  non-homogeneous 
difference  equation.  But  y<  is  a  linear  combination  of  the  first  k  y,’s 
and  of  Cfc+i,. .  and  therefore  Eyiij  is  some  multiple  of  Ee*.  In 
particular,  Eyi€j  will  be  the  same  whether  we  regard  yi,...,yk  as 
fixed  or  as  coming  from  the  stationary  distribution.  This  means  that 
J2|i_y|  is  a  particular  solution  of  the  non-homogeneous  difference 
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equation,  and  so  J2,y  —  R\i-j\  satisfies  the  homogenous  difference 
equation 

Rij  ~  R\i-j\  =  f>ij 

—  ~  H - h  &k6i,j-k)  • 

The  problems  of  summing  the  entire  array  and  of  suming  a  diagonal 
require  different  approaches. 

Summing  the  entire  array  is  very  simple.  Just  use  formula 

(3.1.1)  to  sum  rows  zero  through  Ar  —  1  of  the  array.  Next,  input 

these  sums  as  initial  values  in  (3.1.1)  once  more  to  obtain  the  sum 
of  the  entire  array.  In  our  applications,  the  double  sequences  of 
interest  will  satisfy  the  same  difference  equation  in  both  indexes, 
but  the  above  procedure  works  in  general. 

Next  assume  we  want  to  compute  where  the  dou¬ 

ble  sequence  w  satisfies  the  same  order  k  difference  equation  in  both 
indexes  (t,  j  >  0).  The  difficulty  is  that  the  one  dimensional  se¬ 
quence  (for  fixed  j)  does  not  satisfy  an  order  k  difference 

equation.  We  will  see,  however,  that  an  order  k  vector  difference 

equation  is  satisfied. 

Let  va  =  [ya,a»- ••>!/<», a+fc-i]T  (the  superscript  T  indicates 
transpose),  and  let  F  be  a  k  x  k  matrix  defined  by  the  requirement 
that  f,[ya_fc,...,ya_i]T  =  [ya_*+1,...,ya].  Note  that  to  calculate 
Fv  for  any  vector  v  we  simply  shift  the  components  up,  calculate 
the  bottom  component  from  the  difference  equation,  and  discard 
the  original  top  component.  In  other  words,  we  need  not  store  the 
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entire  matrix  F  in  the  computer,  and  furthermore  multiplication  by 
F  is  much  easier  than  usual  matrix  multiplication. 

Now,  from  the  difference  equation  and  the  definition  of  F  we 
know 

(1)  va  4-  Fva—i  ax  +  •  •  •  +  Fkva-kak  =  0. 

This  is  the  required  vector  difference  equation.  From  this  we  deduce 

OO 

yi  va  —  (I  +  ax  F  +  •  •  •  +  &kFk)  1 

a=0 

x  Vo  4"  (t/j  +  Fvq&i)  H - +  (wfc-j  H - Fk-lvoak-X)  . 

Let  A(z)  =  akzk  H - hao*°-  It  is  easy  to  show  that  the  eigenvalues 

of  F  are  rx , . . . ,  r*  with  eigenvectors  [I, . . . ,  r^-l]r,  where  the  r,  are 
the  roots  of  A(z~l)  =  0,  which  are  all  less  than  one  in  absolute  value. 
(If  the  roots  are  not  distinct  we  can  draw  the  same  conclusions  using 
a  limiting  argument).  Therefore 

k 

I+alF  +  --  +  akFk  =  H(I-riF), 

x 

where  each  factor  is  invertible  since  the  matrix  norm  ||r,-F,||3  is  less 
than  one.  In  other  words,  the  inverse  in  the  expression  for  v<* 
exists. 

In  applications  of  the  above  formula,  we  will  not  generally  know 
vi,  ••  Let  6,  be  defined  as  •  •  >y«,fc-x]r,  for  t  =  0, .. .  ,k- 
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1.  These  6’s  will  be  known  and  because  va  =  Faba,  we  may  rewrite 
(1)  as 

OO 

y!  va  =  {I  +  a\F  4 - +  akFk)~l 

a= 0 

x  6o  +  F(bi  +  toai)  H - 1-  Fk~l(bk~i  H - 1-  6o<ik-i)  . 

The  matrix  polynomials  may  be  efficiently  computed  in  nested  form 
(“Horner’s  rule”),  which  has  the  added  advantage  that  we  only  mul¬ 
tiply  by  F ,  which  is  easy. 

For  convenience  we  will  refer  to  the  function  which  takes  as 
input  y0, . . . ,  y*-i  and  produces  £<T  to  48  output  as 

SUM01(yo,...,yfc_i)> 

and  that  which  produces  y,  as  SUMn(y0, . . . ,  yfc—i) -  The  first 
subscript  (0  or  1)  gives  the  starting  index  for  the  summation,  and 
the  second  indicates  the  univariate  case.  Similarly, 

SUMmI^o, 

0 

and 

OO 

SUMia(6o>  •  •  •  i&fc-i)  =  Vq» 

i 

Difference  recursions  indexed  by  (...,-1,0,1,...) 
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3.4  An  example.  We  want  to  extend  our  arsenal  of  summa¬ 
tion  methods.  The  technique  we  are  about  to  develop  will  greatly 
simplify  certain  calculations.  It  is  best  illustrated  with  a  simple 
example  for  which  it  is  not  really  necessary. 

Suppose  that  a  sequence  {y»  :  t  >  0}  satisfies  the  usual  stable 
autoregressive  difference  equation  of  order  k  and  is  moreover  sta¬ 
tionary  (that  is,  y0,...,yk-i  are  chosen  according  to  the  limiting 
distribution).  We  want  to  calculate  ]C<!-oo  cum(yo,y*)>  which  is 
the  same  as  Yl^oo  &**  where  J2,  is  the  covariance  of  the  process  at 
lag  t.  The  covariances  Ri  do  obey  a  difference  equation,  and  the 
straightforward  method  to  compute  the  sum  is  to  partition  it  into 
two  sums,  £o°  +  Xllio  •  The  latter  sum  is  the  same  as  • 

Though  easy  to  do  in  this  case,  splitting  the  sum  would  be  much 
more  difficult  if  we  were  dealing  with 

OO 

cum(y0yt,y,y,). 

i,y=— oo 

For  one  thing,  we  would  need  to  estimate  k3  fourth  order  cumulants, 
where  k  is  the  order  of  the  autoregression.  This  is  analagous  to  the 
need  to  know  Rq, . . . ,  Rk-i  in  order  to  compute  £3*  Ri.  Further¬ 
more,  although  the  sum  J27js-oo  involves  only  two  indexes,  some 
of  the  component  sums  will  involve  three  indexes  after  shifting.  For 
example,  the  sum  over  (»  >  0 ,j<  0),  after  adding  j  to  each  index 
to  shift  the  indexes  to  non-negative  values  involves  three  indexes. 
This  means  we  would  in  effect  have  to  compute  three  dimensional 
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sums  of  the  form 

00 

i=0 

where  w  obeys  a  difference  equation  in  each  index.  This  all  can  be 
done,  but  it  is  not  the  best  way. 

Let  us  write  a  formal  ^-transform: 

OO 

(1)  <p(zi,zi)=  cum(yj,yi)z{4' 

y,i=-oo 

We  say  “formal”  because  this  sum  does  not  converge  for  any  (zj ,  z3). 
Ignoring  all  such  convergence  problems  for  the  moment,  we  could 
recover  cum(yy,  yj)  from  the  double  contour  integral  (zi  and  za  each 
take  values  on  the  unit  circle): 


1 

(2*t)2 


/ 


<P{Zl,Zl) 

zi+izi2+i 


d(zi,za). 


If  A(z)  =  aitz*  +  •  •  •  +  aQz°,  we  can  “solve”  for  <p  by  multiplying 
both  sides  of  (1)  by  A(zi)A(zt),  where  the  a’s  are  the  coefficients 
of  the  difference  equation: 


<*oyy  •+■••■+  afcy,_fc  =  ey, 

which  we  assume  to  hold  for  all  integer  j,  positive  and  negative. 
Equation  (1)  becomes 

00 

Mzi)A(z7)<p(zltz2)  =  £  cum(€y>ei)zfz£. 

«',j=-oo 


•  •  1.  •  »  » 
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E_2  „m_m 
a(zl  z2  • 

m=—oo 

Substituting  back  into  the  inversion  integral  for  cum(yy,  yi)  converts 
the  double  contour  integral  into  the  sum  of  the  products  of  two 
countour  integrals,  namely 


These  contour  integrals  are  equal  to  t «/_m  and  Wf_m  respectively, 
where  to,-  is  0  for  negative  indexes  and  satisfies  the  difference  equa¬ 
tion  with  coefficients  ao,...,<*fc  subject  to  initial  conditions  w0  = 

l,(u>i  +  aiu/o)  =  •••  =  (u)k-i  H - +  aif-irva)  —  0  (as  follows  by 

writing  the  z-transform  of  this  tt>).  This  suggests 

OO 

cum(yy,y{)  =  ^  u/y_mu;{_m. 


This  conclusion  is  in  fact  valid,  as  we  will  show  later.  Then 


cum(y0,yj)=  ^  w-mwt -mcr\. 


f=  — OO 


/,m=— ao 


Note  that  the  sum  with  respect  to  m  is  a  two  dimensional  re¬ 
cursion,  but  the  sum  with  respect  to  l  is  one  dimensional,  and  is  the 
same  for  all  m.  Summing  first  over  l  and  then  over  m  gives 


OO  /  00  \  3 

y)  cum(yo,yi)  =  )  • 
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Even  the  sum  Wi  is  easy  to  compute.  It  is  1/A(1),  and  so 

oo  2 

H  cum(yo ,yi)  = 

When  this  general  method  is  applied  to  higher  order  quantities, 
the  only  new  parameters  to  be  estimated  will  be  third  and  fourth 
moments  of  e,  as  opposed  to  cum(yo,yt,yy,yj)  for  all  i,  j,  and  l 
between  0  and  fc  —  1.  Furthermore,  we  will  only  have  to  sum  at  most 
two  dimensional  recursions. 

3.5  Validity  of  the  method.  Let  {V,ty  :  (i,j)  G  Z2}  be  a  double 
sequence  on  Z2  which  satisfies  a  difference  equation  in  each  index. 
The  difference  equation  is  not  assumed  to  be  homogeneous.  In  the 
example  of  the  previous  subsection,  cum(y,-,  y y)  was  such  a  Yij.  Let 
B\  be  the  linear  operator  which  shifts  the  first  index,  and  let  be 
that  which  shifts  the  second:  B\  shifts  Yi.y  into  position  (t  +  1,;). 
If  6,,J  is  the  double  sequence  which  is  identically  zero  except  that 
entry  (t,  j)  is  one,  then  B\.6%'3  = 

Multiplying  the  formal  z-transform  of  the  double  sequence  of 
cumulants  {cum(yy,yj)}  by  the  product  AfeijAfa)  to  get  the  for¬ 
mal  z-transform  of  cumulants  {e,tJ}  =  {cum(e,-,  ey)}  corresponds  to 
the  following  linear  operation: 

\Y  =  A(Bl)A(B1)Y  =  e. 

The  linear  operator  A{B\ )  is  by  definition  of  the  polynomial  A  equal 
to  I  +  a\,Bi  +  •  ••  +  In  the  example  of  the  last  subsection, 
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e  =  T°°  We  need  to  show,  in  some  sense,  that  the 

bounded  linear  functional  A  is  invertible  with  a  continuous  inverse. 
Because  A{Bi)A(Bi)  can  be  factored  into  terms  (I  -  r,Bi)  and 
(I  —  r,J5a)  with  jr,j  <  1,  it  will  suffice  to  assume  A  =  /  -  rB,  where 
|r|  <  1  and  B  is  either  B\  or  B2. 

Define  the  measure  u{l}  ~  (1  +  l7)~l  for  l  €  Z,  and  n{l,m)  = 
v{l)u{m).  Consider  Ifin)  on  Z3,  for  p  €  (l,oo).  The  map  A  is  a 
bounded  linear  functional  on  LP{p),  and  furthermore  it  is  one  to  one. 
To  see  this,  observe  that  AT  =  0  and  T  ^  0  together  imply  that  Yjj 
grows  exponentially  fast  as  j  -*  — oo,  so  Y  £  !?((*)•  (This  might 
not  be  the  case  had  we  chosen  ft  to  decay  exponentially  instead  of 
polynomially). 

Next  we  want  to  find  the  inverse  A"1  =  {I~rB)~l.  The  obvious 
candidate  is 

A"1  =  I  +  rB  +  r2B7  H - . 

We  can  bound  the  operator  norm  ||2?fcj|  (=  sup||Bi||/||i||)  by  3k7, 
so  in  fact  the  series  defining  A-1  is  Cauchy  and  the  candidate  inverse 
is  well  defined  and  a  bounded  linear  functional  (hence  continuous). 
Furthermore, 

A~1AT=  lim  (/  +  »•£  +  ■•• +  r*fl,)(/-rfl)T 

<-*oo 

=  lim  (T-r,+lfl,+1T) 

(—*00 

=  Y. 

In  our  applications,  A  will  be  the  composite  of  two  or  more  maps 
of  the  form  A(B),  and  AT  will  have  a  simple  form.  For  example, 
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in  the  last  subsection  A Y  =  e  =  <73<5m,m.  The  inverse  of  6°'° 

is  w  =  {wiWj},  where  u/,-  is  zero  for  negative  indexes,  tt/0  =  1, 

and  (iu!  +  axu;o)  =  •  •  •  =  (tOfc_i  H - +  a^-ituo)  =  0.  By  linearity, 

A ~l6a’0  has  (;',/)  element  Wj-av>i-0-  In  our  example,  Y2-oo 
converges  to  e  in  Lp(/x),  so  continuity  of  A-1  guarantees  that  Yij  = 
72m=-o 0^wi-mwj-m-  (Being  the  Lp(p)  limit  in  this  case  implies 
being  the  pointwise  limit).  Because  the  same  holds  for  any  order 
of  summation  of  X^oo  ,  the  pointwise  convergence  is  abso¬ 

lute.  In  general,  when  Y  is  bounded  A Y  =  e  will  be  bounded,  so 
will  converge  to  e  in  &(/*}.  We  have  proved  the  follow¬ 
ing,  which  clearly  extends  to  any  number  of  indexes. 

3.6  Theorem.  Let  zk-\ - ba*2°  have  a  11  its  roots  in  the  interior 

of  the  unit  circle,  and  let  A(B)  denote  the  operator  (1-1 - +  akBk). 

Let  Y  be  a  bounded  sequence  on  Z3,  and  e  =  A(Bi)A(B-2)Y,  where 
Bi  and  Bi  are  shift  operators  on  the  first  and  second  indexes  as  in 
the  previous  discussion.  Then 

OO 

YlJ  =  Y  €a.0Wl-aWj-0, 

0,0  —  —  OO 


where  w  is  as  in  the  previous  paragraph.  The  convergence  is  abso¬ 
lute.  □ 


Derivation  of 
corrections 


Introduction 

We  will  be  using  Cornish-Fisher  expansions  as  discussed  in  the 
Introduction  to  derive  more  accurate  confidence  intervals  for  the 
autoregressive  process  x  =  {x,  :  t  >  1}  which  satisfies  our  usual 
stable  difference  equation 

(Xi-n) +  ■■■•¥  ak(xj-k  -  M)  =  «/• 

We  assume  the  difference  equation  is  strictly  stable,  which  means 
that  the  roots  of  the  characteristic  polynomial  zk  +  •  •  •  +  akZ°  all 
lie  in  the  interior  of  the  unit  circle.  In  our  derivation,  we  take 
as  given  the  validity  of  the  Cornish-Fisher  expansion  and  of  the 
methods  developed  in  the  preceding  chapters.  Refer  to  the  articles 
by  Taniguchi  (1984),  Abramovitch  and  Singh  (1985),  and 
Gotze  AND  HlPP  (1983)  for  a  further  discussion  of  sufficient  condi¬ 
tions  for  the  validity  of  Cornish-Fisher  and  Edgeworth  Expansions. 
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We  will  assume  that  the  errors  €,■  are  independent  and  identically 
distributed,  have  a  Lebesgue  density  component  which  is  positive 
on  an  interval,  and  have  finite  moments  of  all  orders.  (Recall  that  a 
general  distribution  may  be  expressed  as  the  sum  of  a  distribution 
absolutely  continuous  with  respect  to  Lebesgue  measure,  a  discrete 
distribution,  and  a  singular  distribution.  The  first  of  these  is  the 
Lebesgue  density  component). 

First  we  will  form  a  zero  order  pivot,  which  is  the  usual  pivot 
tn  given  by 

.  _  V*(2  ~  M) 

ln  -  - i  . 

y/Vn 

Here,  we  are  given  £n  is  the  sample  mean,  which  is 

an  estimate  of  the  asymptotic  mean  lim,— «>  Exi  denoted  by  n\ 
and  vn  is  the  estimate  of  the  variance  constant  «  for  2n,  »  = 
limn_*oo  nE(tn  —  n)2.  Our  first  order  correction,  Tj,  is  a  Comish- 
Fisher  expansion  based  on  tn: 

T\  =  fn  +  6n~1/2  +  pt^n~1/2. 

The  Cornish-Fisher  expansion  as  given  in  KENDALL  AND  STUART 
(1977)  is 

Ti  =  tn  —  Iti  —  -K3  (fj  —  1), 

where  «,•  refers  to  the  tth  cumulant  of  tn.  Our  correction  is  the  same, 
except  that  we  must  estimate  the  cumulants.  That  we  estimate 
instead  of  using  the  true  values  does  not  matter  to  the  level  of 
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accuracy  required  for  a  first  order  correction — see  ABRAMOVITCH 
and  Singh  (1985). 

The  next  step  is  to  form  a  Cornish-Fisher  expansion  from  Tj , 
where  now  the  /c,-’s  refer  to  cumulants  of  Ti : 

T,  =  T,  -  K,  -  i(r,2  -  1)  -  i(/c,  -  im 

+  l*i*sr,  -  -3Ti)  +  $(41?  -  7Tt). 

(See  Kendall  and  Stuart  (1977)  or  Hill  and  Davis  (1968) 
for  the  statement  and  derivation  of  the  Cornish  Fisher  expansion). 
T2  differs  from  a  standard  normal  random  variable  by  op(n~l).  It 
will  be  possible  to  obtain  tn  as  a  cubic  polynomial  in  T2,  (which  we 
take  to  be  standard  normal)  instead  of  vice  versa.  This  will  make 
the  formation  of  confidence  intervals  easier. 

The  zero  order  pivot 

4.1  Algorithm.  We  begin  by  specifying  a  possible  pivot  fn,  which 
will  be  used  in  what  follows.  Let  tn  =  n*/3Vn  —  p)>  where  the 
estimate  vn  of  the  variance  constant  for  Sn  is  obtained  as  follows. 

1.  Estimate  the  asymptotic  mean  p  by  p  =  In  =  n-1  i,-. 

2.  Estimate  covariances  Rj  for  ;  =  0, . . . ,  k  by 

1  n~3 

R]  —  —  ~  *»»)(**'+*  —  *n) 

n  •— — 

1=1 

n-j 

J^(x,'  -H  +  H-  in)(xi+j  —  fl  +  p—  In) 

1 


» 
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=  ^  “  M)(*i+S  ~  M)  +  ~  M)5 

71  *  7* 


(M  ~  ®n) 


n 


n-j 


n-j 


L«=l 


«=l 


=  -  53(x*‘  ~  m)(*.-+j  -  m)  +  ~  z*)2 

ft  *-■*  f» 


(*»  ~  M) 


2n(£„  -m)“ 


(1) 


-zo-  53 

n—j+i 


=  ^53(*»-M)(*i+i -/*) 

n  «=i 

-  ^  $3  (*i  -  fi)(xi+j  -  m) 


«=n— y+i 

-  (*.  -  M)3  (i  +  j) 

(*n  ~  H)  '  3 


n 


J3(x,  -/i)+  53  (*<-/*) 

1  n-j'+l 


The  last  form  is  convenient  for  analysis,  while  in  practice  Rj  would 
be  calculated  by  a  method  like  the  fast  Fourier  transform. 


3.  Solve  the  Yule- Walker  equations  for  the  estimated  autore- 
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gressive  coefficients  aT  —  [dt , . . . ,  dfc]T : 


Ro  . 

.  & 
i 

_ i 

i - 

i _ 

Ri‘ 

Rk-i  ■ 

i _ 

-a*. 

. Rk . 

Or,  denoting  the  above  matrices  by  Af,  a,  and  —  R, 

m  =  -it. 

4.  Compute  the  estimate  ~3  of  the  prediction  error  variance 
(variance  of  the  e’s  of  the  mr  ;1): 

&l±<r>  =  lR<)  +  -..  +  akRk 

—  Rq  +  aT  R. 

5.  Estimate  the  variance  constant  v  by 

_ _ _ 

(1 +  dt  +  •  •• +  d*)3’ 

6.  The  pivot  is  tn  =  -  y). 

The  first  order  pivot 

4.2  Introduction.  In  order  to  make  a  first  order  correction 
to  the  confidence  interval  for  the  stationary  mean  y  =  Etx,  we 
need  to  estimate  Et n  and  Et 3,  where  tn  is  the  zero  order  pivot. 
Specifically,  if  Etn  =  crin-1/3  +  a^n-1  +  0(n“3/3)  and  Et*  = 
£ii»“1/3  +  #jn-1  +  0(n-3/3),  the  first  order  pivot  Ti  is  given  by 
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where  9  =  -3c«i/2  +  /3i/6  and  p  =  a i/2  -  /3i/6.  This  follows  from 
the  form  of  the  Cornish-Fisher  expansion  given  at  the  beginning  of 
the  chapter  along  with  the  relationships  between  moments  (/*,)  and 
cumulants  («,): 

*1  =  Mi. 

*  3  =  M3  —  3msMi  +  2m?  • 

In  deriving  T\,  we  will  ignore  terms  which  are  Op(n-1).  (The 
definition  of  “Op(n-1)”  is  similar  to  that  of  “Opfn"1)”:  a  sequence 
{zi  :  i  >  1}  is  said  to  be  Op(n-1)  if  nzn  is  bounded  in  probability). 

4.3  First  moment  of  tn.  From  a  Taylor  expansion, 

t„  =  n~l/2(tn  -  n)  (v'1/2  -  ^u-3/2( vn  -  t/)^  +  Opin'1). 

Then  given  our  regularity  conditions  stated  in  the  first  paragraph  of 
this  chapter  (i  =  {i,- :  *  >  1}  satisfies  a  stable  order  k  autoregressive 
difference  equation  whose  errors  e,-  are  independent  and  identically 
distributed  with  positive  Lebesgue  density  component  on  an  interval 
and  with  moments  of  all  orders), 

Etn  =  En'1l2{tn-n)v'1l 2 

-  !*!-«/»*-*/»(*  -  v)(*„  -  M)  +  Oin'1). 

The  regularity  conditions  imply  that  E{sn  —  m)  =  n-1M-1(^*)  + 
0(n-3/3),  where  as  usual  (i~1{  A*)  means  53“  E(x{  —p).  Therefore 
the  first  element  of  the  n~ 1!2  coefficient  of  Etn  (that  is,  aj )  is  simply 
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fi~l  (A2n)v-1/2.  Our  analysis  is  conditional  on  *i, . . . , (k  is  the 
order  of  the  difference  equation)  and  thus  from  the  formulas  for 
summing  difference  recursions, 
v~l,2ii~l( A2)  =  trl/3SUM0i(ri  ifc  -/i) 

-i/a  (»i  ~  A*)!1  +  •  •  •  +  Qfc-i)  H - -I-  (gfc  -/>)■! 

(!+  •  +  <**) 


and  the  natural  estimate  is 

-\H  (*1  _  *n)(l  H - 1-  Ofc-i)  H - !■(**“  *n)  •  1 


The  other  element  of  ari  is 

tr3/Vi(A* At/)  i  En(tn  -  M)(-l/2)v-3/3(vn  -  v). 

Expanding  vft  —  w  as  a  Taylor  series  in  (a2  —  a2)  and  £V(d,-  —  a<) 
gives  the  first  order  approximation 


(1)  vn-v  = 


a3 -a3 
(1  H - h  afc)3 

_  2<r2 

(H - (-a*) 


y(di  —  Oj  H - 1-  a*  —  Ofc). 


The  autoregressive  coefficients  d  satisfy  Md  =  —  R,  while  for  the 
true  values  we  have  Ma  —  -R.  Writing 

(M  +  AM)(a  +  Aa)  =  -(R  +  A  R) 

shows  that  AMa  +  MAa  +  AMAa  =  -  Ail  or 

a  A 

a  —  a  =  Aa 


=  -(M  +  AM)"1  (Ail  +  AMo) 
=  -M"1  (AR  +  AMo). 
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We  may  similarly  find  A (M  x)  =  M-1  —  M  1  by  writing 
(M  +  AM)  (A/”1  +  A(M“1))  =  I 

=>  A(M“l)  =  -M-'AMAT1 
=  -M-1  AMM-1 


to  first  order.  If  desired,  then,  we  have  a  second  order  approximation 
for  Aa: 

(2)  Aa  =  -M~l(AR  +  AMa) 

+  M“‘  AMM“*  (AR  +  AMa). 


For  the  present,  we  only  need  the  first  term  on  the  right  hand  side. 

The  remaining  part  of  the  n"1/3  coefficient  of  Etn  (that  is,  ari), 
is  then 


2 


fi)v-3/'{vn-v) 


■r  V3/2n 


E(xn-n)(a'-cr*) 


(!  +  •••  +  a*)3 


2<73j E(xn  -  /*)lr(a  -  a) 

(1  +  •  •  •  +  afc)3 

_  >  -3/a-[g(*n  ~°7) 

2  "I  (l+...  +  afc)3 

_  2<r2E(Xn  -  /t)lT(-l)Af~l(AR  +  AMa)' 
(1  +  •••  +  a*)3 


Recalling  that  the  notation  fi~l(A2AR)  denotes  the  most  signifi¬ 
cant  vector  coefficient  of  E(tn  —  n)AR,  with  similar  definitions  for 
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^-‘(AiAJZo)  and  the  matrix  ^"‘(AiAAf),  linearity  shows  that 

the  formula  above  is 

1  _3/J E[tn  -  n){c7  -  a7)n 
2V  (14 - 1-  a*)a 

+ - i — - (f*"1  (Ax  A/2)  +  M-H^xAMifl)]. 

(14 - h  a*)3  J 

Next  we  must  express  a 2  —  a7.  We  have 
&7  =  R0  +  arR 

=  (Ro  +  A  Ro)  4-  (a  +  Aa)r(/2  4-  A  R), 

so  to  first  order 

A<t3  4  a*  -  a7 

=  Ai2o  +  A  a?  R  +  a3"  A/2 
=  A/2o  -  [A/'1  (A/2  4-  AAfo)]T/2  4-  aT A/2. 

Because  a  =  —M~ 1 R  =  — M~T R ,  this  is 

A(73  =  A/2o  +  (A/2t  4-  orAM)a  +  aT  A/2 
=  A/2o  4-  2orA/2  4-  arAMa. 

Therefore  to  first  order 

^“‘(AxA^3))  =^_1(AiA/2o) 

4-  2a3V-l(A£A/2) 

-(-  aT fi~l(AiAM)a. 

Combining  results  shows 
"y/nv~3f7ji~l(^  Av)  = 
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’/i-1(A*AJ2o)  +  2  aT  n~l(AxAR)  +  ar/i-1(AiAAf)a 

.  ^ 


2a3ITM-1(M_1(AiA£)  +/i"1(AzAM)a) 

+  ^(i)3 


Therefore,  the  n~1/3  coefficient  of  Etn  is 

-1/3  (si  _  *n)(l  H - b  flfc-l)  4 - 1"  (Xfc  -  Xn)  •  1 

(3)  0,1  =  v~  - <i 

1  _3/2  fi~l(Ax&Ro)  +  2aT fi~l(AxAR)  +  aT ii~*  (AxAM)a, 

~2V  y/"[  A(l)2 

2a3lTM-1(/*_1(AxAi2)  + /i-1(AxAAf)a) 

+  3(7)3 


We  are  left  to  compute  n~l(AxA  Rj)  for  j  =  0,...  ,k.  Recall 
the  form  of  Rj  given  in  point  2,  equation  (4.1.1),  in  the  section  on 
the  zero  order  pivot. 

Let  R'j  =  n~l  £“(x,-  “  M)(*i+y  ~  /*)•  Then 
M~1(axA  R'j'j  =  hJ1(AxA  R^ 

oo 

=  Yi  E»(xo  -#»)(*<-/*)• 

*=— OO 

The  results  of  Chapter  2  show  that  n~x(AxA Rj)  =  /i-1  ^AxAjR'). 
It  is  worth  noting  in  passing,  however,  that  the  end  effects  in  the 
computation  of  Rj  can  be  significant,  though  they  are  not  in  this 
instance.  In  particular,  the  results  of  Chapter  2  also  show  that 
fi~7(Ax2ARj)  ?£  n~7(^Ax7ARj'j. 
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Toestimate  J2-oo  E*(xo-p)(xj-p)(xi-p),  we  use  the  methods 
of  Chapter  3.  This  illustrates  the  utility  of  Theorem  3.6.  We  have 

E,(xi  ~  li)(x j  -  n)(ii  -  n)  =  cum(xj  -  /i,xy  -  /*,  x,  -  n) 

OO 

=  53  cum(e,e,€)ym_jym_yym_< 

m=— oo 


where  the  cumulants  are  with  respect  to  the  stationary  distribution 

and  y  has  z-transform  l/.4(z)  =  (a0z°  H - h  a^z*)-1.  This  means 

that  y,-  is  zero  for  negative  indexes  and  obeys  the  autoregressive 

difference  equation  with  yo  =  1,0  =  yi  +  aiyo  =  •  •  •  =  yk-i  H - h 

a*_,y0.  Therefore 


«=  —  00 


53  Et(xo  -  fi)(Xj  -  fi)  ( i(  -  n)  = 

00  oo 

K3(«)53y*  H  VO-mVi-i 

«=0  rr»=  — oo 
_  «3(c)  V' 

“  ■JTjj’  Zv  y mym+j- 

'  '  m=  0 


We  already  know  how  to  compute  [£  ymym+y]y_o  >  because  this 
is  just  SUMo2[yay^]'  The  special  structure  of  \yaye>\  facilitates  the 
calculation  of  the  sum,  as  reflected  in  the  algorithmic  summary  at 
the  end  of  the  chapter. 

4.4  Third  moment  of  tn.  For  the  first  order  pivot  Ti  we  also 
need  to  estimate  J3lt  the  n~lf7  coefficient  of  Et*,  or 
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=  n3/2(2n  -  /i)3  ^«"3/2  -  ^-5/2(v„  -  +  OpCn"1). 

Therefore  /?!  =  v~3/2ji"2( Ax3)  —  3ti-5/,2/i-2(Ax3Av)/2.  Now 

OO 

li~7( Ai3)  =  £.(*o  -  /*)(x,  -  /i)(iy  -  m) 


|,J  =  — OO 


—  ^  1  *3^)  )  ]  y-myi-mt/j-r 


i,j=  — oo  m=  — oo 


_  *a(«) 

From  the  fact  that  /x-2( Ax3)  =  /u72(Ax3)  +3^~1(Ai)«i  we  find 

_  K3(«)V-3/2  „ 

Pi  =  -~1)a--+3M-,(Ai)ti-1/a 

—  ^w-5/2^~2(Ax3Aw) 

it 

K3(e)v"3/2  _1/3 

=  A( ip— +  3^  l(Ax)V  »/* 

—  ^v_3^2/*-l(AxAv). 

it 

Replacing  all  quantities  in  the  above  equation  and  in  equation  (4.3.3) 

for  a!  by  their  natural  estimates  enables  us  to  compute 

*  3  .  1  - 

»  =  -j«i  +  gA, 

-  1  - 

'  =  2°l  “  6* 

T  -  i  4.  *  . 

i  1  —  tfl  -r  -7=  +  — =. 

v71 


The  second  order  pivot  T2  derives  from  a  Cornish  Fisher  ex¬ 
pansion  of  Ti .  The  expansion  is 

t2  =  r,  -  icx  -  l-K3(T?  - 1)  -  ±(k7  - m 
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+  -  ^<C4(^’l3  ~3Tl) 

+  ^*?(U't3-77'1)+Op(n-1), 


where  is  the  *th  cumulant  of  Ti.  It  is  not  hard  to  show  that 
and  k3  are  o(n-1),  while  k2  and  *4  are  0(n-1).  Because  tn  —  Tx  is 
we  will  rewrite  the  above  as 


T2  =  Ti  —  ~(k2  -  1)<„ 

-  -3t„)  -I-Opln-1). 


Again,  the  fact  that  we  estimate  n«a  and  rwc4  does  not  matter  to 
°p(n-1)— -see  Abramovitch  and  Singh  (1985). 

4.S  Second  cmnulant  of  T\ .  We  have 


Kt-  1  =  ET?  -  E^Ti  -  1 
=  ET?  -  1  +  o{n~l) 


=  4"+^  +  w)’-,+0,n~,, 

* i  +  3£  +  *S. 

n  n  n 
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+  2pP'  +  ^ 

n  n 

+ - Vn - + - ^ - +  “(n  >' 

Estimation  of  most  of  these  quantities,  and  of  those  needed  for  /c4, 
is  straightforward.  The  necessary  calculations  are  summarized  in 
the  subsection  at  the  end  of  the  chapter,  along  with  a  few  brief 
explanations  and  references  to  Chapters  2  and  3.  In  the  current 
subsection,  we  will  discuss  points  which  require  further  elaboration, 
and  will  point  out  some  pitfalls. 

In  many  cases  we  will  need  the  most  significant  coefficient  of 
covariance  between  the  point  estimate  xn  and  the  product  of  two 
estimates,  pq.  In  other  words,  we  want  the  most  significant  part  of 
E(£n-n)(pq-pq).  To  “first  order”  the  change  pq-pq  is  pAq+qAp, 
so  we  expect  the  coefficient  to  be 

PH~l(AxAq)  +  qn~l(AxAp). 

This  conclusion  is  generally  correct,  but  one  must  be  a  little  careful 
of  the  logic.  The  true  change  includes  a  A pAq  term,  but  the  ex¬ 
pectation  E(AxApAq)  is  of  lower  order  than  the  other  terms.  If, 
however,  we  were  considering  E(tn  —  p)2(pq  —  pq),  we  would  need  to 
account  for  the  extra  EAx2ApAq  term,  because  this  is  of  the  same 
order  as  EAx2Ap  and  EAx2Aq.  In  short,  one  must  carefully  apply 
the  results  of  Chapter  2. 

We  now  turn  to  specific  items,  following  the  same  order  as  in 
the  second  order  summary  at  the  end  of  the  chapter. 
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I.  /»72 ^Ax2  .  From  the  definition  of  iZ'  and  the  mixed 
third  moment  identities,  this  is 

OO 

Y  cum((i0  -  fi){xj  -  n),  ( xi  -  /i)(xm  -  n)). 

l,m=  —  oo 

Again  we  will  see  the  usefulness  of  Theorem  3.6. 

Let  =  cum((x,  -  (i){x3-  -  n),  (n  -  fi){xm  -  fi)).  Apply¬ 

ing  the  autoregressive  difference  equation  to  each  coordinate  trans¬ 
forms  this  sequence  to  the  following  e  sequence: 

=  cum(€,ey,€tem) 

(  Ee*  —  a*  for  »  =  j  =  l  =  m, 

<r4  for  *  =  l  ^  j  =  m, 

I  a*  for  »  =  m  ^  =  1, 

'  0  otherwise. 

Theorem  3.6  then  shows 

OO 

—  (Ee*  —  3tr4)  }  '  y«-piO'-pyi— pl/m— p 

p—  —  oo 

00 

+  O*  Y  yi-vVi-^yi-pym-q 

p,c/=  — oo 
oo 

+  <?*  Y  yi-pyj-iyi-iym-p. 

p,<j=-oo 

Recall  that  y  is  the  sequence  whose  z-transform  is  1/A(1).  Including 
those  pairs  (p,q)  for  which  p  =  q  in  the  latter  two  sums  is  compen¬ 
sated  by  the  -3<t4  in  the  first.  Summing  the  above  over  (l,  m)  G  Z3 
now  leads  to  the  result  stated  in  the  summary. 
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2.  /*72(A2 n&Rj)-  The  expression  of  the  summa--y  is  obvious, 
except  for  the  -jpJl(&x*ARj)  term.  This  comes  from 

1  n 

-E(£n  -  m)2  *  —  (*i  -  /*)(*.+ 3  -  /*)• 

n  n— 7+1 

To  find  the  value  of  this,  it  is  only  necessary  to  equate  0(n?)  terms 
of 

0{n)  =  cum(£"(x,-  ~  aO>5Zi  [x*  -  A»)>  (xn-i  ~  /*)»  (*n-i+j  ~  p))- 

3.  p~l  (aRj-ARI ^ .  Again  use  Theorem  3.6.  From  the  moment 
identities,  we  want  to  estimate 

oo 

cum((x0  -  -  p),  (xm  -  p)(xl+m  -  p)) 

m=— oo 

under  the  stationary  distribution.  Theorem  3.6  reduces  part  of  this 
to  a  sum 

oo 

y-pyj-pym-pyt+m-p- 

m,p ——oo 

This  accounts  for  the  S}Si  term  appearing  in  the  summary.  The 
rest  is 

oo 

y-py:-qym-Pyi+m-q 

m,p,q=—oo 

00 

+  2^  y-pyj-qym-qyi+m-p- 

m,p,q=-ao 
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We  will  illustrate  the  method  for  the  first  of  these  two  sums, 
sum  over  p  and  then  over  q.  The  result  is 

OO 

m=— oo 

For  convenience  suppose  l  >  j.  Then  both  and  are 

equal  to 

OO 

m= 0 

In  general,  this  leads  to 

0<m<|I— y| 

except  that  if  l  =  j  we  need  to  subtract  S3  from  this  (because  then 
Er  and  overlap). 

4.  /i-2(Ai3Aa).  Note  tha.t  we  consider  the  entire  expres¬ 
sion  for  Aa  given  in  the  first  order  correction,  Formula  (4.3.2). 
In  general, we  will  need  to  consider  higher  order  errors  to  compute 
£-3(A23A())  or  1  (A(-) )  for  any  “(•).” 

5.  /i-1(AxA(63)).  The  estimate  of  Ee*  is  given  by 

Et*  =  n~l  ^  cf 
fc+i 

n 

=  n~l  ^2  [(*<  -*«)  +  •••  +  h{Xi-k  -  *n)]3  • 
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This  reduces  to  our  problem  to  estimating  the  most  significant  co¬ 
efficient  of 


and 


E{i»-p)n  1  XX 


E(in  -  n)n  1  XX?  -  c?). 


The  first  of  these  coefficients  we  refer  to  as  /*7l(A*nA(n_1 
Because  e3  —  e3  is  a  polynomial,  we  could  write  the  second  exactly 
as  a  Taylor  series,  but  we  only  need  the  following  approximation: 

ft 

£(*n-M)n-1  £(€?-«?)  = 


3 n~l  XI  «?  — A{1)(*„  -  /*)  +  X)  Aoi(I«-J  ~  M) 
1  3=1 


Part  of  this  is,  to  our  order  of  aproximation, 

-ZA(\)c7E(tn-n)2. 


The  rest  is 

k  n 

X;  EAaj(£n  -  A»)3 n~l^2€i(xi-j  ~  /*)• 
i- i  »=i 

The  inner  sum  estimates  cum(e3,z,_y)  ,  which  is  zero.  The  expec¬ 
tation  EAa}  (£n  -  fi)A  (cum(c3,z,_y))  is  of  lower  order,  though  of 
course  the  estimate  need  not  be  zero. 
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6.  Calculation  of  M3  and  /^(Az^AM?).  The  algorithm  of 
the  summary  uses  Horner’s  rule  or  nested  evaluation  of  the  matrix 
polynomial  defining  M3.  From  the  stated  algorithm,  it’s  clear  that 
M3  is  computed  as 

F  (•  •  •  (FaicI  +  ik-il)  H - )  4-  aol. 

The  calculation  of  £7l(A*nAM3)  parallels  this. 

The  rest  of  the  calculation  of  k,7(Ti)  is  straightforward. 

4.6  Fourth  cumulant  of  7\ . 

1.  /i~3  ( AlJ).  From  Theorem  3.6,  the  fourth  order  cumulant  of 
the  sum  (x*~  A*)  under  the  stationary  distribution  is  #c4(c)/A( l)4 . 
Now  equate  the  O(n)  terms  of 


«4  =  M4  -  3M33, 

where  the  k  and  p’s  are  moments  and  cumulants  of  (*<  —  A*) 
under  the  stationary  distribution. 

2.  /*“ 1  ( AiZy).  By  this  we  mean  the  n_l  coefficient  of  E(Rj  — 
Rj).  Note  that  in  general  in  computing  £_2(A())  we  must  account 
for  higher  order  terms.  For  example,  to  compute  /t_I(Aa)  we 
use  the  first  and  second  order  terms  of  Ao  given  in  the  first  order 
correction. 

3.  The  expression  given  derives  from  the  relationship 
between  non-central  moments  and  cumulants, 


«4  =  A»4  ~  4^3A*i  -  W  +  12M3M13  ~  6/ii4, 
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which  in  the  case  of  J\  reduces  to 

ETf  -3-6(K2(ri)  -  1). 


Inverting  the  correction 

We  now  discuss  how  to  invert  the  first  or  second  order  correc¬ 
tion.  Suppose 


9fl  Pn^n 

Z  =  tn  +  -~+  ^  +  -=-2-, 

y/n  y/n  n  n 


and  define 


g{x,0,p,v,u)  =  - 


By  9n  (and  so  on)  we  mean  the  estimate  of  9  based  on  *i , . . . ,  x„. 
Tlien  g{tn,9n,pn,vn,un)  is  Op(n~l),  as  are  g(£,  9n,  pn,vn,  un)  and 
tn  —  £.  Ignoring  terms  which  are  Op(n~3/7)  and  writing  g  for  the 
function  <?(£,  9n,  pn,  vn,  wn)  and  /  for  the  derivative  of  g  with  respect 
to  its  first  argument, 

0  V  +  0„(n-3/3) 

=  J  +  <Ms  +  (i»-«)^)  +  0p(n-3/3) 

=  g  +  9,g  +  op{n-3/2) 

\  y/n  y/n  n  n  / 

+  ^  (29./^  +  2^«3)  +  Opln-3/’). 
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Therefore 


tn  —  £  ~  r-  (On  +  />n£3) 

Vn 

+  -  (-t'ni  +  2 9npnt  -  un£3  +  2 pi?)  +  Op{n~3l3). 
n 

Here  (  is  regarded  as  a  normal  random  variable.  This  equation 
enables  us  to  convert  normal  quantiles  (£  values)  into  tn  quantiles 
directly,  rather  than  solving  for  the  tn  quantiles. 

Algorithmic  summary 

In  this  section,  we  summarize  the  computations  necessary  to 
make  the  first  and  second  order  corrections.  As  usual,  xt-  is  an 
autoregressive  process  of  order  k.  The  autoregressive  coefficients  are 
oo  =  1,  , . . . ,  at,  the  errors  are  e*, . . . ,  en,  the  variance  of  the  e’s  is 

a3,  while  v  is  the  asymptotic  variance  constant  for  Xn.  The  analysis 
is  conditional  on  Xj,. ..  ,  x*.  The  first  and  second  order  corrections 
depend  on  the  method  of  estimating  v,  and  we  assume  that  the 
method  detailed  at  the  beginning  of  the  chapter  is  used. 

4.7  First  order  summary. 

1.  Compute 

M  =  SUMoi(li  —  Sfit  •  ■  •  >  Xk  —  t  n). 

2.  Estimate  e,-,  i  =  k  +  1, . . . ,  n.  The  estimates  are 

=  (*i  —  *n)  +  —  *») 

•  "  +  Ofc(x,_fc  —  xn). 
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3.  Calculate  j&e3  =  n_1  J2k+i  • 

4.  Calculate  1/j4(1)  =  (1 H - ha*)-1.  This  is  the  sum  l/«» 

where  y,-  has  the  2-transform  1/4(2).  Recall  that  A{z)  =  0*2*  -I- 
- ha0z°  (o0  =  1). 

5.  Calculate  0 y«y*'+j>  (3  =  0, ...,k  —  1)  as  follows.  First 
solve  the  triangular  system 


A 

=  «!  = 


a0J  |y*-i 


Recall  that  F  is  a  shift  matrix  defined  by  the  requirement  that 


The  matrix  F,  then,  involves  the  estimated  autoregressive  coeffi¬ 
cients.  The  desired  sums  are  the  solution  of 


(l+  -  +  akFk) 


.  S  yalfa+k—l 


6.  Compute  £  yayq+  *  from  the  sums  of  the  previous  paragraph 
using  the  fact  that  the  sums  £a>0  yaya+»  as  a  function  of »'  obey 
the  same  order  k  difference  equation. 

7.  Compute  the  covariances  between  the  point  estimate  tn  and 
estimated  covariances  Rf. 

1  ( yalfa+jt 

a>0 
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for  j  =  0, . . . ,  k. 

8.  Let  ii'1  (A£nAM)  denote  the  symmetric  k  x  k  Toeplitz 
matrix  with  (t,j)  entry  equal  to  A  7 1  (  A2„  A/t|,-_  y| ) .  Similarly 
let  fiJi(A±nAR)  be  the  column  vector  with  tth  entry  equal  to 
fiJi(A2nARi),  for  t  between  1  and  k. 

9.  Compute  the  estimated  coefficient  of  covariance  between  the 
point  estimate  I»  and  the  vector  of  estimated  autoregressive  coeffi¬ 
cients  a: 


A.  ‘(Al#Aa)  =  -M  1  (fcl(AxnAR)  +  ^(A^AM)*)  , 

and  also 

A7l(A*„A(lra))  =  lrAr‘(A£nAa). 

10.  Calculate  the  estimated  coefficient  of  covariance  between 
the  point  estimate  and  prediction  error  variance  estimate: 

A.-'fAj.Af*1))  =  Ar‘(Al„Afl„)  +  ATfi7‘(  Ai„Aa) 

+  aTAr‘(A*„Afl). 

11.  Calculate  the  estimated  coefficient  of  covariance  between 
the  point  estimate  and  vn  : 


fi~l(AlnAv) 


ji;'(A£nA(<T')) 
(!  +  •••  + a*)3 


2(72At  ^A^A^a)) 

(1  +  --  +  C*)3 


12.  The  estimate  of  on,  that  is,  the  n-1/3  coefficient  of  Etny  is 
given  by 


di  =  A  ‘(AS*)**1/3  -  ^Ar1(A*nAu)v“3/3. 


TTTT 
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13.  The  most  significant  coefficient  of  the  third  central  moment 
of  the  point  estimate,  under  the  stationary  distribution  is  then: 


.4(1)* 


14.  The  estimated  n-1/2  coefficient  of  Et*  is  then 

-  |Ar‘(Ai„A«)v73/J. 

15.  One  may  now  compute  the  coefficients  for  the  Cornish- 
Fisher  polynomial: 

-  3-  U 

*”~2a!+6/Jl’ 

.  1.  U 

^=2ai'6^* 

16.  The  first  order  corrected  pivot  is  then 

Ti  =  tn  +  9n~xl2  +  ptln~1/7, 

and  one  may  compute  the  p-quantile  of  the  distribution  of  tn  by 
zp  -  9n~ll2  -  pz7n~ll7, 

where  zp  is  the  p-quantile  of  the  standard  normal  distribution. 

4.8  Second  order  summary.  The  bulk  of  this  calculation 
consists  of  the  calculation  of  the  second  cumulant  of  T%.  Though  it 
may  appear  lengthy,  the  calculation  is  not  too  computer-intensive. 
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Calculate  /i“3(Ax3)  using  the  moment  identities: 

OO 

i 

2.  A,'3  (A*;)  =  -2 .  SUM.,  (ET  ft.  • .  ■ .  s?  ft)- 

Estimate  the  difference  between  the  above  and  its  non-stationary 
equivalent.  From  Chapter  2,  this  amounts  to  summing  the  double 
array  6  below  over  the  entire  first  quadrant.  ‘ 

3.  Si,j  =  (x,-  —  £n)(*y  —  *n)  — 

for  1  <  i,j  <  k. 

4.  TEMP,-  =  SUMoi($,-,i, . . . , 

for  i  = 

5.  A-J(A*J)-A,-3(Ai;)  = 

SUM01(TEMP,,. . .  ,TEMPt). 


6.  A-'(A4)  =  Ar3(^i3)  +  (A-’(AiJ)  -  ■ 


In  the  next  sequence  of  equations  we  will  compute  £-3(Ai3  A-Ry). 
Recall  that  R'  is  n~l  £){x,-  -  /t)(x,+y  -  p). 


7. 


• _ 2  {  a  a  m  \ 


iTe4  —  3<74  r-' 
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Refer  to  the  first  order  correction  for  ^2yaya+jt  which  has  already 
been  computed. 

9.  (Tl  (^i2y)  =  SUMoa  [(x*  -  2n)(xi  -  2n)  -  £|»-i|]y » 

where  the  value  for  j  =  k  is  obtained  from  the  values  for  j  = 
0, . . . ,  k  —  1,  using  the  fact  that  as  a  function  of  j,  these  quantities 
obey  the  autoregressive  difference  equation.  Next,  for  j  =  0, . . . ,  k, 
compute 

10.  A7«(A*»A*,)  =  -3  (A7‘(AiJ))’  . 

-  + a.-j(a*;a  *j), 

u.  A-2(Af;Ajj,)  =  A7i(Ai;Aj?,) 

+  2p-‘(  AiJAr1  (Ai.Aii;.) 

In  the  following  formulas,  5,-  will  stand  for  ]Ca>o  SfaVa+it  which 
has  been  computed  in  the  first  order  correction.  Let  S3i|/_y|  stand 
for  the  sum  £a>0  S<*Sa+ |»-y|-  Compute: 

12.  ‘S'a.ji-yi  =  SUM02(5’a5^)|(_J|, 

and  iterate  the  above  out  to  S-j^k  using  the  autoregressive  differ¬ 
ence  equation  (which  5a,,-  obeys  as  a  function  of  *).  Next  estimate 
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+  a 4  j^252i|j_;|  -f  25ati+j 

+  5Z  s*S\t-j\-i 

0<«'<|i-j| 

+  ^*Si+i-i  • 

0  <i<t+j 

The  estimate  of  n~x  ^A R'ARj^  is  then 

.  .  (  TEMP  -  S$crA,  if  /  =  i  #  0 

14.  A 71  f  A5J  A R\)  =  l  TEMP  -  2Sq&4,  if  /  =  ;  =  0 

'  [  TEMP,  otherwise. 

15.  fcl{&RjARi)=jxJ1(&R'JAR'iy 

In  order  to  estimate  p~3(A£*  Aa) ,  first  compute  the  following: 

16.  A72(A£’AJ2,  A#)  = 

A71  (&Rj&Ri ) 

+  2H-1(AtnARj)fiJl(A£nARt). 

Let  be  a  k  x  k  matrix  whose  (*,  l )  entry  is  one  if  |»  —  l\  =  j 
and  zero  otherwise;  for  i  <  j  <  k  let  RW  be  a  column  vector  with 
elements  1, . . . ,  k  which  are  zero  except  for  the  jth  element  which  is 
one;  and  let  R(°)  =  0.  Compute: 

17.  Wjl)  =  M~l(Rl})  +  M(j)a), 

18.  wjV 

J 

k 

K'iAaAR,)  =  -  T,w}"ji7'(ARlAIli), 


19. 
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for  /  =  0, . . . ,  k. 

20.  A_2(Ax*Aa)  = 

o<»',i<fc 

21.  A72( AilAaiARi )  s  ji-l(Atl)jiJl{AaiARi) 

+  2a7'(A2nAai)d7'(AtnARi), 

22.  ii~7{AxlA(a2))  =  ^{AilARo)  +  aT  (Ai'AR) 

+  RTiT7  (Ai*  Aa) 

k 

+ 53  ft*  (Ai» A/z») » 


23. 


24. 


ft 

A7l(AaAa,)  = 

1=0 

Jig 

A7l(^(lro)3)  =  A7l(Aa,Aay), 


i,J=l 


25.  A7l(^3)^dra))=E 


1=1 


A7l(Aa,A.R0) 


k 


H-E^MAa.Aa,) 

j=i 

+  53ajA7l(Aa,Ai?y)l 
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26.  A73(A**A(<r2)A(lra))  = 

d7'{A£2)fc'{A(a2)A(l  Ta)) 

+  2  fc1  (AxnA(lT  a))^1  (AinA(<r2)) 

27.  ^73^Af2A(lra)2)  = /i71(A23)Ar1(^(ira)3) 

+  2(A71(A*«A(lra)))2 

28.  ,-Wv)  m 

23 *fi-*(A£lA(lTa)) 

Mi)3 

_  2Ar2(^^nA(g2)A(ira)) 

Mi)3 

K’UtlWTa)’) 

A(  l)4 

Next  we  will  compute  (i72  {Ai^Ay2) . 

k 

29.  j*-l(A(<;2f)  =  [aia3fiJl{ARiARj) 

»'.j=o 

+  RiRjiiJl(AaiAaj) 

+  2a,J?//i7I(AayA.R,)), 


where  A7l(A<*oA())  —  0  because  a0  =  1  is  known. 

^‘(a^2)3)  4d4A71(^(iTa)3) 

30.  *71(Av2)  =  - \ - L  + - ^ - L 

^  K  1  A(l)*  il(l)« 
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ftjl(A£nAM7) «-  A71  (A£nAa*.)  J 

£71(AxflAa  o)  =° 

For  l  —  k  —  \  down  to  0, 

l»r1(A*fiAMa)4-/i71(A*(tAf,)Ma  +  #,A7l(A*.AA#a) 
+A71(A2flAa<)/ 

A/2  < —  FMi  +  cl(I. 

36.  «-‘(A*„A(A/f‘))  =  -M2lK‘(Al^M2)M^' 


M*  = 


ak-l  •••  ®0 


The  coefficient  of  covariance  of  the  point  estimate  with  Af3,  namely 
ftJl(A£nAMz),  is 

’A  71(Ai„Aa0)  0 

38.  : 

A  71(A#rtAofc_l)  ...  A71(A2nAo0), 

The  diagonal  elements  in  the  matrix  above  are  all  zero,  but  they  are 
written  as  they  are  to  indicate  the  structure  of  the  matrix. 

39.  A7l(A*nA  (M3-1))  =  -A/3-lA71(AJ„AM3)M3-1 

40.  5  =  [50,...,5fc_llT 

=  [Erf.-.E»»-<-‘] 

41.  A71(AinA5)  =  A7l(A**A 

+  M3-iA71(AJnA(M3-1))ei. 


•  .*  >  V  *.»  V  *  *  V  V  V  *A-  V.S  V  •  *»  '  •  k  »  *»  *  •  »  *.  r. 
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_  4<72£71(A(<r2)A(lra)) 

i(i)5 

31.  /i72(Ai2  At/2)  =  A71  (Ai2)A71(^v3) 

+  2(A71(^2»Av))3. 


This  enables  us  to  compute 

ra(Aita)  _  A-a(A^A») 
“  «Ar‘(A*i)  n^r'CAii))’ 

Ara(A^Ana) 

»(A.-‘(Aia))3' 


33.  Ar‘ (Ai„A(n-‘  ££?))  =  SUMo,  (o, . . .  ,0, £<*), 

34.  A71(a*„A(&»))  =A,-‘(Ai„A(n-‘  £€,)) 

-3#aA(l)Ar‘(A*i). 

The  next  series  of  computations  will  give  us  Mj  =  {I  +  aiF  + 

- hdfcFfc)  and  The  matrix  fiJ1(&xn&F)  below 

is  a  kx  k  matrix  identically  zero  except  in  the  last  row: 


£7‘(A*nAF) 


0 

.-£7l(A*BAofc) 


0 
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In  the  following  calculations  we  will  compute  coefficients  of  co- 
variances  between  the  point  estimate  and  estimated  coefficients  from 
the  first  order  correction. 

42.  = 

5Ar‘(^.A(£«4)) 

i(l) 

i(l)» 

.  Ar‘(Ai!.A S)6* 

id) 

43.  Ar,(A*»AMr‘(AinASt))  = 

k  ( 

Y  [  -A71  (^*n  Aoi  JAr1  ( A2„  A-Rfc-/ ) 

1=1  ' 

-  aiAr1  (AfnA(A*7l (A2n A/fc_/))) . 

Define  A71  (A*BA/i7l (A2„A/2))  in  the  obvious  way,  that  is  by 
[A7 1  ( Ain  A/17 1  ( Ain  AJZ, ) )  ] *_  j ,  and  similarly  define  the  coefficient 
A71(A2nA/i7t(A2nAAf)). 

44.  A7l(A*„A(Af-1))  =  -M-1A7l(A2nAAf)M~l. 

45.  A7l(A*nA#i7‘(A*nAa))  = 

-Ar‘(A*.AM-‘)k-'(Al.Afl) 
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+  /*7l(A*nAA/)a 


A7l(A2nA^i  '(&£nAR)) 


+  H7'(Asn&n7'(A£nAM))a 
+Ar,(Ai.AM)Ar'(A*.A»)l 


46.  A7l(A*„AMr‘(A*.A((TJ)))  = 

/i;‘(Al„A,i7'(Al1,Afi0)) 

+  Ar'(A*.Afl1,)Ar‘(A*„Aa) 

+  RT  iiJl(A£nAn7l(^£nAa)) 

+  A71  (A*n  Aar)A7l  (A*n  Afl) 

+  dTAr‘  (A2»  A^7‘  (A**, AJ2)) , 

47.  /*71(AinA^71(A*nAu))  = 

A71(Aa?nAM7  l(AgwA(<ra))) 

MD2 

_  4A71(A2nA(0-3))A71(A2nA(lro)) 
6*a  (A71(AgftA(lra)))2 

i(l)4 

2<73A7l(A2«AM7l(A*nA(lra))) 
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wT  = 


[1  +  dj  +  ■  ■  •  +  dfc_i, . . . ,  1  +  dj,  l| 


49.  ft,  1  (AinAwT)  =  [1  +  dj  H - +  ,  1  +  dj,  lj 

(-l)A.-‘(Af,A(lTa)) 


l 

•=l 


50.  1(A£n))  = 

-SUM0l(l,...,l)Arl(Art) 

+  /tf1  (AinAtl>r)  •  (*1  “  *«i  •  •  •  I  *fc  -  £nlr- 

51.  fi-*(AfHAat)  =  A71(A*»A*r1(A*j)i£,l/a 

-±v;3/*M-l(A£n)M7,(A*nA'') 

~  \v~3f2iiJl(A£n/\n-l(Mn^v)) 

+  ^n5/3  (/“rH^n^W)3* 

£rl(A2nA(£€3)) 

52.  A71  (A..A*7’(A*i))  =  Z 

^VrlCAg^AU^a)) 

i(l)4 
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53.  )  =  /t71(A2nA/i73(AiJ))t>-3/3 

-  |^6/aM7a(A*i)#ir1(A*. Av) 

+  SAr1  (Ain  Am~1  (A*  J)  w“; l'a 
“  £<£3/aA-1  ( A^JAr1  ( Ai„ Av) 

-  ^Ar^^^nAMrMAfnAv))^-3/3 
+  7^5/2(Ar1(A3:„Av))a. 


54.  Ar‘(AJnA0)  =  -^Ar1(A*nAot1)  +  ^Ar!(AiCnAA) 

55.  A71  ( A2nAp)  =  ^Ar1  ( A*rt  Aorx )  -  gA7l  (A*. A/?x ) 

56.  A71/2(Atn3A/>)  =3A7l(A2nAp)v"l/3 

57.  A71/2(Af*A0)  =  Ar‘(A2nA«)V71/3. 


This  enables  us  at  last  to  compute  the  estimated  second  order 
cumulant  of  2\,  namely  k2{Ti): 

2A.-1/3(Atn3A/0 


58. 


ki(Ti)  —  1  —  Et\  —  1  + 


n 


,  2A7I/3(AtwAg)  |  fla  t  3pa 
n  n  n 

+  2gdx  +  2ffix  2gp 
n  n  n 
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Now  we  turn  to  the  fourth  order  cumulant  of  Ti . 

59.  k^(e)  =  Ee4  -  3<74 

60.  £73(Af4)  =  4^  +  6£7‘(A5£)£r2(Aj’) 

61.  a-^ajj)  =  £r3(A*i) + 4a-‘(a*„)A73(a*;) 

+6a.-‘(a*;)  (r3(Ajj)  -  a73(a*;)) 

By  £-1(A Rj)  we  mean  limfl_00  n(Rj  -  Rj).  For  /  =  0 
compute 

62.  A'1  (AS,)  =  A-1  (as;)  -  ,S,  -  A7‘  (AlJ). 

63.  £-‘(Aa)  =  -M-1  (£-l(Ai2)  +£“1(AAf)a) 

k 

+  £  A7‘(ASiAS,)tvg), 

«\j=0 

64.  /x-^ACa2))  =  £-1(Ai2o)  +  ar£~l(A.ft) 

k 

+  i2r£-1(Aa)  +  ^£r1(  AoiAUk), 


£_1(Av) 


£-x(A(<r2)) 

i(l)2 


2<y2Ir£~l(Aa) 

i(iP 


65. 
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66.  £“3(A2jAt;)  =  4jiJ2  (<i£3n)  {&inAv) 

+  6£-2(A22  Av)v* 
-3u2£-l(Aw), 

67.  £73(A*<A«,2)  =3(£7HA*3n))2A7l(Av3) 

+  12£7l(A*n)  (£7l(A2„Av))3, 
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i\£H£ 
nL  vl 


2£-3(Ar*At;) 


3£73(Ag*Av2) 


69.  £-1/2(A/n3A^)  =£72(A23A^)v73/2 

=  3£71(A2rtAtf)V71/2. 

In  the  following,  note  that  £-l/2(Afn3)  is  the  same  (by  defi¬ 
nition)  as  4i>  and  £-l/2(Af»)  is  the  same  as  dj. 

70.  £-l/2(A^8)  =  10£-l/2(A«„3)  -  I5£-1/2(A<«), 

71.  £-l/2(Afn5Ap)  =  I5£7l/2(Azn  Ap)t/7l^a, 

72.  £°(Atna)  =  15. 
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73.  ET*  —  3  =  Et\  —  3  +  ^  40£-1/2(A*„3) 

+  4£-1/3(Afn3A0)  +  4pA~1/3(Afn5) 

+  4/i-l/3(Afn5Ap)  +  603  +  90p3  +  36 9p  . 

And  therefore  the  fourth  order  cumulant  estimate  for  T\  is 

74.  *4(ri)  =  (ET*  -  3)  -  6(«a(Ti)  -  1). 


The  quantiles  tp  for  tn  may  be  estimated  from  corresponding  stan 
dard  normal  quantiles  zp  by 


A  (1  —  p)  confidence  interval  for  n  =  Et  x  is  then 
77.  (2n  -  fl-i/3„»/3f|_p/3(  *n  +  n-i/awi /2tp/a). 
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4.9  A  note  on  validation  of  algebra.  If  for  a  given  autore¬ 
gressive  model  one  inputs  the  true  covariances  and  moments  of  the 
residuals  into  the  algorithms  for  the  first  and  second  order  pivots, 
one  can  find  the  true  values  of  all  the  quantities  estimated  in  this 
chapter,  provided  the  calculations  here  are  correct.  Even  for  the 
zero  order  pivot,  using  true  covariances  will  yield  the  true  asymp¬ 
totic  variance  constant  v.  This  was  done  for  an  AR3  model,  and 
the  theoretical  results  thus  obtained  were  compared  with  simulated 
values.  Simulated  values  are  of  course  the  correct  ones,  except  for 
the  error  of  estimation.  In  this  experiment,  there  were  1,000  data 
points  per  replication  and  10,000  replications.  Here  are  a  few  results 
from  the  second  order  correction. 


Quantity 

Theoretical 

Simulated 

/i-2(Ax2A.Ro) 

22.23  •  101 

22.35  •  101 

H~'(AR 33) 

83.62 

81.06 

Hjl(Aa3AR3) 

-20.45  •  10"1 

-19.76-10- 

/x“2(Ax2Aoj) 

-25.94 

-23.13 

Hjl  (Ax An,1  (ArA.Ro )) 

-28.65 

-29.76 

H~l(AxAiiJl(AxA{<x7))) 

-24.0 

-25.02 

pJ3(AxAAv2) 

32.56  •  104 

32.86  •  104 

(i~3(  Ax*Av) 

36.83  •  103 

48.90  •  103 

ET i4  -  3 

15.90-  10"3 

-7.14-  10“ 

The  above  figures  are  typical.  Note  that  in  the  last  two  rows 
above,  the  agreement  may  not  be  as  good  as  one  might  expect,  but 
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the  numbers  are  in  fact  acceptable.  ET *  -  3  is  one  of  several  sec¬ 
ond  order  quantities  with  an  extremely  high  coefficient  of  variation. 
Because  of  this,  it  is  not  possible  to  obtain  a  good  estimate  given 
current  computing  constraints.  Estimates  of  odd  order  moments  like 
/i”3(Ax4Av)  also  tend  to  have  a  high  coefficient  of  variation.  In 
summary,  almost  all  of  the  estimates  obtained  confirm  the  algebra 
of  this  chapter,  and  none  of  them  arouses  suspicion. 


Numerical  results 


Introduction 

The  usual  test  statistic  tn  is  asymptotically  normal,  as  are  the 
corrected  statistics  T\  and  7?.  The  distributions  of  the  latter  two 
statistics,  however,  converge  more  quickly  to  the  standard  normal 
distribution.  (See  for  example  ABRAMOVITCH  AND  SINGH  (1985) 
and  the  references  cited  there). 

To  understand  the  data  to  follow,  it  may  be  helpful  to  give  an 
informal  review  of  how  we  form  the  usual,  “zero  order”  confidence 
interval  based  on  the  standard  test  statistic  tn,  and  how  we  form 
the  first  and  second  order  corrected  confidence  intervals. 

The  statistic  tn  converges  weakly  to  a  standard  normal  random 
variable.  If  zp  denotes  the  p-quantile  of  the  standard  normal,  the 
approximation 

P[z. os  <  tn  <  z. 95}  «  .9 

becomes  more  accurate  as  n  increases  (“«”  means  approximately 
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equal).  In  fact,  the  error  is  0(n“1/3).  The  event  that  the  true 
asymptotic  mean  /i  lies  in  the  nominal  90  percent  confidence  interval 
constructed  from  tn  is  precisely  the  event 

{■8.05  <  tn  <  z, 95}, 

while,  for  example,  the  event  that  n  lies  above  the  upper  90  percent 
confidence  interval  bound  is  the  same  as  the  event 

0»  <  s.os}* 

This  is  evident  by  rewriting  tn  as  >/n(2n  -  n)/arn. 

Ti  is  a  quadratic  polynomial  in  tn,  and  is  a  cubic  polynomial 
in  tn ,  which  we  might  indicate  by  writing  Ti(tn)  and  For 

values  of  tn  of  interest  and  for  sufficiently  large  n,  these  polynomials 
are  approximately  the  identity  function: 

rl(t„)  =  tn  +  op(n-1/2) 

T2(tn)  =  Tl(tr,)  +  Op(n~l). 

These  polynomial  transformations  may  be  regarded  as  transforming 
the  distribution  of  tn  into  a  (more)  normal  distribution.  Equiva¬ 
lently,  these  transformations  map  quantiles  of  the  distribution  of  tn 
into  corresponding  approximate  normal  quantiles.  Not  that  tn,  7\ , 
or  Tj  are  actually  normal,  but  each  of  the  following  statements  is 
more  accurate  than  its  predecessor: 

^{■8.05  <  tn  <  z. 95}  »  .9, 

P{*.05  <  Ti  <  z_ 95}  «  .9, 

P{z, 05  <  Ta  <  z. 95}  »  .9. 
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In  fact,  the  errors  are  0(n-1/3),  o(n-1/3),  and  o(n-1),  respectively. 

The  inverted  transformations  may  be  written  as  T^l(zp)  and 
T;l{zp).  Instead  of  converting  fn-quantiles  to  approximate  normal 
quantiles,  we  are  doing  the  reverse:  converting  normal  quantiles  to 
approximate  in-quantiles.  The  inverses  of  the  quadratic  Ti()  and 
of  the  cubic  T?(-)  are  not  polynomials,  but  we  have  shown  how 
to  approximate  these  inverses  with  a  polynomial  up  to  the  desired 
order  of  accuracy.  These  inverted  polynomials  are  then  T and 
They,  too,  are  nearly  the  identity  for  fixed  values  of  their 
arguments  in  the  range  of  interest: 

rf'w  =  z+op("~1/J) 

rr'w  =  rr'w  +  o,(n-'). 

We  could  base  modified  confidence  intervals  on  the  statements 

P{z.os  <  Ti  <  z95}  «  -9 
P {^.05  <  Tj  <  ^.95}  w  -9, 

but  conversion  of  these  statements  into  statements  about  (in  —  p) 
requires  solving  a  quadratic  or  cubic  equation.  It  is  therefore  more 
convenient  to  base  intervals  on 

P{Trl(z.os)  <tn<  ‘(z.m)}  «  .9 

P{T3-l(z05)  <  fn  <  Tfl(z.95)}  w  .9. 

The  event  that  the  true  mean  p  lies  above  the  lower  90  percent 
confidence  interval  bound  obtained  from  the  first  order  corrected 
interval,  for  example,  is  the  same  as  the  event  that  tn  <  TT^z.os). 
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In  this  chapter  we  will  compare  the  distribution  of  the  usual 
test  statistic  tn  with  the  standard  normal  distribution  and  with  the 
estimated  distribution  of  tn  obtained  from  the  inverted  versions  of 
7\  and  Tj,  which  we  denote  by  Tf 1  and  Tf1.  As  we  have  said,  in 
practice  one  does  not  use  Ti  or  Ta  directly  to  form  a  confidence  inter¬ 
val  because  this  would  require  solving  a  cubic  or  quadratic  equation. 
Instead,  we  use  the  inverted  corrections  discussed  above  and  given  in 
equations  (4.7.16)  or  (4.8.76).  However,  the  comparison  of  T\  or 
with  the  standard  normal  is  qualitatively  similar  to  the  comparison 
of  the  inverted  statistics  with  tn. 

It  would  be  possible  for  us  just  to  report  the  true  coverage  prob¬ 
ability  for  various  intervals.  In  other  words,  with  what  frequency 
does  the  true  asymptotic  mean  fall  above  a  90  percent  confidence 
interval,  or  below  it,  or  in  the  upper  half,  or  in  the  lower  half?  There 
is,  we  feel,  a  more  informative  way  to  compare  the  various  pivots 
using  what  we  will  call  p-p  plots  which  are  a  variation  on  q-q  plots. 

Suppose  we  wish  to  compare  two  distributions,  Fi  and  Fa.  A 
q-q  plot  plots  the  two  quantiles  corresponding  to  a  given  probability. 
Thus  the  plot  includes  points  of  the  form  (Ff1  (p),  F^fp)).  If  Ft 
is  standard  normal,  and  Fa  is  an  empirical  distribution,  this  is  well 
known  as  a  useful  way  to  test  normality.  It  has  the  virtue  of  testing 
the  tails  (which  are  typically  of  most  interest  to  us)  and  also  the 
virtue  that  if  Fa  is  non-standard  normal  the  plot  will  still  be  a 
straight  line. 

For  our  purposes,  the  p-p  plot  is  better.  Here  we  plot  points 
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(Fi(q),  Fj(g)).  Both  axes  of  these  plots  consist  of  the  interval  [0, 1]. 
Suppose  Fi  is  standard  normal  and  F\  is  the  distribution  of  tn,  and 
that  the  point  (pt ,  p 3)  appears  on  the  graph.  This  means  that  the  p\ 
quantile  of  tn  is  the  same  number  as  the  pa  quantile  of  the  standard 
normal,  which  in  turn  means  that  the  true  probability  that  tn  is  less 
than  the  pa-normal  quantile  is  p\.  The  true  coverage  probability  for 
a  90  percent  confidence  interval  based  on  tn  is  then  pugh  —  Plow 
where  the  points  (phigh,0.95)  and  (piOw>0.05)  appear  on  the  plot. 

Referring  to  the  first  graph  for  Model  1,  four  increments  in 
normal  probability  subdivide  [0,  l]  on  the  y-axis  (“Normal  p”  axis), 
namely  0.05,  0.45,  0.45,  and  0.05.  The  corresponding  increments 
are  indicated  along  the  fn  p  axis,  which  in  this  case  are  0.18,  0.47, 
0.32  and  0.03.  This  means,  for  example,  that  the  probability  that 
tn  is  less  than  or  equal  to  the  0.05  normal  quantile  is  actually  0.18. 
Furthermore,  the  event  that  tn  is  less  than  the  0.05  normal  quantile 
is  the  same  as  the  event  that  the  true  mean  /*  lies  above  the  upper 
90  percent  confidence  limit.  For  reference,  the  dots  on  the  graph 
indicate  the  points  (.1,  .1),  ...,(1,1). 

The  results  of  testing  the  inverted  corrections  are  displayed  in 
a  similar  way.  The  inverted  first  order  statistic,  Tfl,  gives  for  each 
replication  an  estimate  of  any  quantile  of  tn.  This  estimate  of  the  pth 
quantile  is  T1_l(zp),  and  it  changes  with  each  replication.  Referring 
to  the  second  row  of  the  tables  for  Model  1,  the  tn  p  coordinate 
corresponding  to  the  0.95  T”1  p  coordinate  gives  the  frequency  with 
which  tn  is  less  than  or  equal  to  its  estimated  0.95  quantile  obtained 
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from  the  first  order  correction,  namely  T^l(z.3$).  Again,  the  event 
that  tn  is  less  than  or  equal  to  this  estimate  is  the  same  as  the  event 
that  the  true  mean  p  lies  above  the  lower  90  percent  confidence 
interval  bound  as  given  by  the  first  order  corrected  method.  In  this 
way,  we  see  that  for  n  =  200  data  points,  the  probabilities  that 
a  nominal  90  percent  confidence  interval  covers  the  true  mean  are 
0.79,  0.81,  and  0.87  for  the  zero,  first,  and  second  order  methods, 
respectively. 

Models  tested 

We  will  test  the  following  four  models. 

(1)  Autoregressive  model  (r<  -  10)  -  0.5(x,_i  —  10)  —  0.3(x,_a  — 
10)  -  0.  l(x,'_3  -  10)  =  e,  -  1,  where  the  e,’s  are  exponential 
with  mean  1.  The  three  initial  values,  X|,  xa>  and  x3,  are  set 
to  seven,  which  is  three  less  than  the  true  mean.  Results  are 
shown  for  n  =  200  and  n  =  400  data  points. 

(2)  The  same  model  as  (1),  but  with  geometric  residuals:  F{e,-  = 

/}  =  for  j  >  0. 

(3)  The  waiting  time  process  in  the  M/M/1  queue,  {W,-  :  *  >  1}, 
with  traffic  intensity  p  =  0.5.  We  model  this  as  an  AR5  process 
and  set  the  initial  value  to  0.  Results  are  shown  for  n  =  1,000 
data  points. 

(4)  A  Markov  chain  on  the  nonnegative  integers  with  transition 
probabilities  Pi,i+i  =  1  -  Pi,i-i  =  1/3  for  t  >  1.  This  is  a 
discrete  analog  of  the  queue  length  process  of  the  third  model, 
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which  we  also  model  as  AR5.  The  initial  value  is  set  to  zero. 

Results  are  shown  for  n  =  1,000  data  points. 

Note  that  the  third  model  violates  the  independence  assump¬ 
tion  on  the  prediction  errors  and  also  the  assumption  of  finite  au¬ 
toregressive  order.  We  therefore  would  not  necessarily  expect  good 
results  for  this  model.  In  the  second  and  fourth  models,  either  the 
prediction  errors  or  the  process  itself  has  a  lattice  distribution,  which 
again  means  we  would  not  necessarily  expect  good  results. 

Data 

We  have  chosen  to  do  10,000  replications,  a  number  somewhat 
larger  than  those  used  in  JOHNSON  (1978),  GLYNN  (1982a),  and 
Jow  (1982).  Suppose  we  want  to  estimate  the  probability  that  tn 
is  less  than  or  equal  to  the  0.05  quantile  of  the  standard  normal 
distribution,  z0 s-  Our  experiment  is  one  of  binomial  trials  with  p 
approximately  equal  to  0.05.  A  95  percent  confidence  interval  for 
the  required  probability  will  then  have  halfwidth  of  0.0043,  or  about 
half  a  percent.  This  is,  we  feel,  a  desirable  level  of  accuracy  for  such 
an  experiment.  The  halfwidths  corresponding  to  z.so  and  z0i  are 
about  1.0  percent  and  0.2  percent. 
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5.5  Discussion  of  data.  The  graphs  indicate  that  the  method 
performs  very  well  for  both  of  the  true  autoregressive  models  tested. 
This  is  to  be  expected,  especially  in  the  case  of  the  first  model  which 
satisfies  all  of  the  sufficient  criteria  posed  in  Chapter  4.  In  the  lat¬ 
ter  two  models,  we  also  see  some  improvement  despite  the  depen¬ 
dence  of  residuals,  but  we  would  expect  to  see  more  improvement 
if  this  dependence  were  taken  into  account.  Using  the  Fast  Fourier 
Transform,  it  is  possible  to  perform  the  first  order  correction  in  the 
case  of  dependent  residuals  in  0(n  log(n))  time,  but  a  second  order 
correction  would  be  more  difficult.  This  is  one  possible  area  of  fu¬ 
ture  research.  Autoregressive  models  have  the  virtue  of  tractability 
and  of  being  a  convenient  way  to  take  into  account  some  kind  of 
nonstationary  behavior.  But  in  many  cases,  some  other  method  of 
modeling  this  nonstationarity  may  be  more  desirable,  and  the  basic 
ideas  and  tools  developed  here  can  be  used  to  obtain  asymptotically 
more  accurate  confidence  intervals  in  these  other  contexts. 
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